{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robust nonlinear regression in scipy\n",
    "TAGS: Optimization and fitting\n",
    "\n",
    "AUTHORS: Nikolay Mayorov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the main applications of nonlinear least squares is nonlinear regression or curve fitting. That is by given pairs $\\left\\{ (t_i, y_i) \\: i = 1, \\ldots, n \\right\\}$ estimate parameters $\\mathbf{x}$ defining a nonlinear function $\\varphi(t; \\mathbf{x})$, assuming the model:\n",
    "\\begin{equation}\n",
    "y_i = \\varphi(t_i; \\mathbf{x}) + \\epsilon_i\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where $\\epsilon_i$ is the measurement (observation) errors. In the least-squares estimation we search $\\mathbf{x}$ as the solution of the following optimization problem:\n",
    "\\begin{equation}\n",
    "\\frac{1}{2} \\sum_{i=1}^n (\\varphi(t_i; \\mathbf{x}) - y_i)^2 \\rightarrow \\min_\\mathbf{x}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Such formulation is intuitive and convinient from mathematical point of view. From the probabilistic point of view the least-squares solution is known to be the [maximum likelihood](https://en.wikipedia.org/wiki/Maximum_likelihood) estimate, provided that all $\\epsilon_i$ are independent and normally distributed random variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So theoretically it is not optimal when $\\epsilon_i$ have distribution other than normal. Although in engineering practice it is usally not important, i.e. if errors behave as some reasonable random variables with zero mean a result of least-squares estimation will be satisfactory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The real problems start when data is contaminated by outliers (completly wrong measurements). In this case the least-squares solution can become significantly biased to avoid very high residuals on outliers. To qualitatively explain why it is happening, let's consider the simplest least-squares problem:\n",
    "\\begin{align}\n",
    "\\frac{1}{2} \\sum_{i=1}^n (a - x_i)^2 \\rightarrow \\min_a\n",
    "\\end{align}\n",
    "And the solution is the mean value:\n",
    "\\begin{align}\n",
    "a^* = \\frac{1}{n} \\sum_{i=1}^n x_i\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now imagine: $a = 1, n=4, x_1 = 0.9, x_2 = 1.05, x_3=0.95, x_4=10 \\text{ (outlier) }$. The solution $a^* = 3.225$, i.e. completely ruined by a single outlier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the well known robust estimators is l1-estimator, in which the sum of absolute values of the residuals is minimized. For demonstration, again consider the simplest problem:\n",
    "\\begin{equation}\n",
    "\\sum_{i=1}^n |a - x_i| \\rightarrow \\min_a\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the well-known solution is the [median](https://en.wikipedia.org/wiki/Median) of $\\{x_i\\}$, such that the small number of outliers won't affect the solution. In our toy problem we have $a^* = 1$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only disadvantage of l1-estimator is that arising optimization problem is hard, as the function is nondifferentiable everywhere, which is particularly troublesome for efficient nonlinear optimization. It means that we are better to stay with differentiable problems, but somehow incorporate robustness in estimation. To accomplish this we introduce a sublinear function $\\rho(z)$ (i.e. its growth should be slower than linear) and formulate a new least-squares-like optimization problem:\n",
    "\\begin{equation}\n",
    "\\frac{1}{2} \\sum_{i=1}^n \\rho \\left((\\varphi(t_i; x) - y_i)^2 \\right) \\rightarrow \\min_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Turns out that this problem can be reduced to standard nonlinear least squares by modifying a vector of residuals and Jacobian matrix on each iteration, such that computed gradient and Hessian approximation match the ones of the objective function. Refer to the [paper](http://link.springer.com/chapter/10.1007%2F3-540-44480-7_21) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implemented choices of $\\rho$ are the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Linear functior which gives a standard least squares: $\\rho(z) = z$.\n",
    "2. [Huber loss](https://en.wikipedia.org/wiki/Huber_loss): $\\rho(z) = \\begin{cases} \n",
    "z & z \\leq 1 \\\\\n",
    "\\sqrt{z} - 1 & z > 1\n",
    "\\end{cases}$\n",
    "3. Smooth approximation to absolute value loss, \"soft l1 loss\": $\\rho(z) = 2 (\\sqrt{1 + z} - 1)$\n",
    "4. Cauchy loss: $\\rho(z) = \\ln(1 + z)$.\n",
    "5. Loss by arctan: $\\rho(z) = \\arctan z$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions 2 and 3 are relatively mild and give approximately absolute value loss for large residuals. The last two functions are strongly sublinear and give significant attenuation for outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loss functions above are written with the assumption that the soft threshold between inliners and outliers is equal to 1.0. To generalize, we introduce the scaling parameter $C$ and evaluate the loss as \n",
    "\\begin{equation}\n",
    "\\hat{\\rho}(r^2) = C^2 \\rho \\left( \\left(\\frac{r}{C} \\right)^2 \\right)\n",
    "\\end{equation}\n",
    "\n",
    "Note that if residuals are small, we have $\\hat{\\rho}(r^2) \\approx \\rho(r^2) \\approx r^2$ for any $\\rho(z)$ defined above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate the behaviour of the introduced functions we build a plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from matplotlib import rcParams\n",
    "rcParams['figure.figsize'] = (10, 6)\n",
    "rcParams['legend.fontsize'] = 16\n",
    "rcParams['axes.labelsize'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r = np.linspace(0, 5, 100)\n",
    "\n",
    "linear = r**2\n",
    "\n",
    "huber = r**2\n",
    "huber[huber > 1] = 2 * r[huber > 1] - 1\n",
    "\n",
    "soft_l1 = 2 * (np.sqrt(1 + r**2) - 1)\n",
    "\n",
    "cauchy = np.log1p(r**2)\n",
    "\n",
    "arctan = np.arctan(r**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1058c1ef0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAGECAYAAABptmcuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VNXaxuHfSgKhhRYQBSSAfBQVROk9IiJ2sQA2FMtB\nkWNFRT3SVGwHUVEsKAoWjh0sIKgYCCQIgqigoFISBOm9BJLM+v7YmWQmJJA22TOZ576uuTLZs2fm\nHUflYa39rmWstYiIiIhIcIpwuwARERERyZ/CmoiIiEgQU1gTERERCWIKayIiIiJBTGFNREREJIgp\nrImIiIgEsaAKa8aY+saYucaYlcaYX40x/846PtIY87cxZlnWrY/btYqIiIiUBhNM66wZY04ETrTW\nLjfGVAGWApcC/YF91trnXC1QREREpJRFuV2AL2vtZmBz1v39xpjfgXpZDxvXChMRERFxSVBNg/oy\nxjQEWgM/ZB0aaoxZbox5wxhTzbXCREREREpRUIa1rCnQj4G7rLX7gYlAY2tta5yRN02HioiISFgI\nqmvWAIwxUcCXwCxr7Qt5PB4HfGGtbZXHY8H1YURERESOwVp73Mu8gnFkbTLwm29Qy2o88LocWJHf\nk621uoXobeTIka7XoJu+v3C86bsL7Zu+v9C73X+/pUePgo8vBVWDgTGmC3At8Ksx5ifAAg8D1xhj\nWgMeYD0w2LUiRURERIrorbfgs89g0SKoVatgzwmqsGatXQhE5vHQ16Vdi4iIiEhJmj8fHnzQ+Rkb\nW/DnBeM0qISp+Ph4t0uQYtD3F7r03YU2fX+hYc0a6NcP3nsPmjcv3HODrsGgOIwxtix9HhEREQl9\ne/ZAp04wdCgMGZJz3BiDLUCDgcKaiIiISIBkZMBFF0GTJvDSS/6PFTSsaRpUREREJEDuvhusheef\nL/prBFWDQaA1bNiQlJQUt8uQEhAXF8f69evdLkNERCRfL70E338PSUkQVYzEFVbToFnDjaVYkQSK\nvksREQlms2bBTTc5Qa1Ro7zPKeg0aFiNrImIiIgE2ooVcMMNznpq+QW1wtA1ayIiIiIlZMsWuPhi\nGD8eunQpmddUWBMREREpAWlpcNllcP31cO21Jfe6umZNQpK+SxERCSYejxPQrIVp08Ac90o0XbMm\nIiIiUmpGjoT162Hu3IIFtcLQNGiIGzVqFBEROV9jREQEY8aMcbEiERGR8DJlirON1IwZULFiyb++\nRtZCnDEG4xPhFy1aRP369V2sSEREJHwkJMD99zs/TzghMO+hsFbGtG/f3u0S8nXkyBHKly/vdhki\nIiIlYvVq6N8f3n8fTj01cO+jadAyJvc0qHea9K+//uKiiy4iJiaGhg0b8thjjx313O3bt3PbbbdR\nv359KlSoQIsWLZg0aVKe5zRr1ozKlSvToEEDrr32WjZt2uR3nvd9V65cSZ8+fYiJiaF///6B+dAi\nIiKlbPt2uPBCeOIJ6NUrsO+lkbUyzjtFevnllzNo0CDuvfdevvjiC0aOHEmDBg244YYbANi3bx9d\nunTh8OHDjBkzhoYNGzJ79mxuv/12jhw5wh133AHAzp07iY6OZuzYsZxwwgn8888/jBs3jq5du7Jq\n1arskTPv+1522WXcfPPNDB8+3O/aOhERkVB1+DD07QtXXAG33BL491NYCwPGGIYNG8bAgQMB6Nmz\nJ9999x3Tpk3LDmvPP/88GzZsYMWKFTRu3Dj7vF27djF69Ghuv/12IiIiaNq0KS+88EL2a3s8Hjp3\n7kyDBg2YNWsWl156qd/73nXXXQwdOrQUP62IiEjgeDwwaBDUqQNPPlk676mhjmMwJnC30nbBBRf4\n/X766aeTmpqa/fvs2bPp0KEDcXFxZGZmZt969+7Njh07+O2337LPfeWVV2jdujUxMTFERUXRoEED\njDGsXr36qPe97LLLAvehREREStmIEbBuHbzzDpTWhJFG1o6hLK25WrNmTb/fo6OjSUtLy/5969at\nrFmzhnLlyh31XGMMO3bsAGDChAncddddDBs2jN69e1OjRg08Hg8dOnTwez2vk046qYQ/iYiIiDsm\nT3YWvE1ODswSHflRWBMAYmNjqVOnDi+++GKeOwM0a9YMgA8++IBevXrxzDPPZD+2fv36fF/XuDGM\nKCIiUsK+/RYeegjmzQvcEh35UVgTAPr06cNLL73EySefTK1atfI97+DBg1SrVs3v2OTJkxXKRESk\nzFq5Eq65Bj76CJo3L/33V1gTAO655x4+/PBDunbtyj333EOzZs04cOAAq1atIjExkenTpwNOqHvm\nmWd48sknad++PXPnzuXjjz92uXoREZHA2LzZWaLjueegRw93alBYKwN8R7Vy72iQ+/H8jletWpWk\npCTGjBnDM888w8aNG6levTrNmjXjiiuuyD5vxIgR7Nmzh+eff560tDTi4+OZM2cOjRs3LvD7ioiI\nhIIDB+Dii+Gmm+C669yrw+R1fVKoMsbYY32erN3tS7EiCRR9lyIiEkiZmc5aarGxTmNBIMYfsv4s\nO+4ra+kOERERER/Wwp13wqFD8Prr7iy55UvToCIiIiI+xo2DxETnlseKVqVOYU1EREQky0cfwQsv\nOGup5Vr8wDUKayIiIiLAwoVwxx3wzTdQv77b1eTQNWsiIiIS9v74w9mY/Z134Iwz3K7Gn8KaiIiI\nhLUtW+D88+GJJ+C889yu5mgKayIiIhK2DhyAiy6C66+Hm292u5q8aZ01CUn6LkVEpLgyMuCyy5y9\nPt98s/SX6NA6ayIiIiL5sBaGDHEC22uvub+W2rGoG1RERETCztix8OOPMG9ecKyldiwaWQtxo0aN\nIiIiAo/HUyKvFxERwYgRI0rktURERILRlCkwaRJ89RXExLhdzfFpZC3E5bVxu4iIiORtzhx44AFI\nSICTTnK7moLRyJqUuvT0dLdLEBGRMLRsGVx7LXzyCbRo4XY1BaewVkasXbuWiy66iJiYGBo2bMhj\njz2W/djbb79NREQEqampfs/xTqHmZq1l7NixnHzyyVSqVIkePXrw888/H3Xep59+SqdOnahcuTI1\natSgX79+bNiwwe+cRo0acf311/PWW2/RokULoqOjmTlzZgl9ahERkYJZu9ZZouO116BrV7erKRyF\ntTLAWsvll1/OOeecw4wZM+jbty8jR45kypQpQP5TpfkdnzJlCrNmzeLll19mypQpbNmyhV69erF7\n9+7sc1599VWuvPJKTj/9dD755BNef/11VqxYQXx8PAcOHPB7ve+//57x48czatQovv76a1q1alXC\n/wRERETyt3079OkDjzwCl1/udjWFp2vWygBjDMOGDWPgwIEA9OzZk++++45p06Zxww03FPr10tLS\n+Oabb6hQoQIA7du35//+7/8YP348o0eP5sCBAwwfPpybb76ZSZMmZT+vffv2NG3alDfffJM777wz\n+/ju3bv56aefqF27djE/qYiISOF4F7298kpn389QpLB2DGZ04C7ctyNLdkHXCy64wO/3008/neXL\nlxf5tbxBDSAuLo6OHTuSnJwMQFJSEvv27eOaa64hMzMz+7x69erRvHlz5s+f7xfWOnbsqKAmIiKl\nLiMDBgyAZs2craRClcLaMZR0oAqkmjVr+v0eHR1NWlpakV6rTp06eR777bffANi2bRvWWs4555yj\nzjPGHFXLSaHSbiMiImWGtXD77ZCeDm+8EdyL3h6PwloY8I6SHTlyxO/4jh078jx/y5YteR6rV68e\nALGxsQBMnTqVU0899ahzY3ItWqOlRUREpLSNHAk//QTffx/8i94ej8JaGIiLi8Nay4oVK2jSpAkA\nmZmZzJkzJ8/zZ86cyaFDh6hYsSIA69evZ9GiRTz88MMAdO7cmZiYGP7880+uu+660vkQIiIiBfTK\nKzBtGixcGBqL3h6PwloYaNeuHaeccgr3338/mZmZREdHM3HiRA4fPpzn+RUrVqR3794MGzaMtLQ0\nRo4cSfXq1bn77rsBZ+Ts2WefZejQoWzdupXzzz+fatWqsXHjRubNm8fZZ5/NgAEDSvMjioiIAPDp\np/D445CY6GzQXhYorJUB+U0zeo9HRkby+eefc8cddzBo0CBq1qzJ3XffTYcOHRgzZsxRz7nhhhuo\nVKkSQ4cOZceOHbRv356PPvqI6tWrZ5/3r3/9iwYNGvDss88ybdo0MjIyqFevHt26daN169Z+r6dp\nUBERKQ3z5sFtt8Hs2dC4sdvVlBxjbehcRH88xhh7rM9jjKEsfd5wpu9SRER8/fornHOOM/2ZR/9b\nUMr6s+y4IxpaFFdERERCWkoKXHABvPhi6AS1wlBYExERkZC1bRv07g3DhjlrqpVFCmsiIiISkvbv\nhwsvdHYnuOsut6sJHF2zJiFJ36WISHg7csTZRiouDl5/PTQXvS3oNWsKaxKS9F2KiIQvjweuvRbS\n0uCjjyAqRNe2KGhYC9GPJyIiIuHIWrj7bti0yVmiI1SDWmGEwUcUERGRsmLsWJg/31lTLWs3xTJP\nYU1ERERCwmuvweTJsGABVKvmdjWlR2FNREREgt5HH8GYMc6o2kknuV1N6VJYExERkaA2Zw4MHQrf\nfAOnnOJ2NaVP66wJAJMnT6Zp06ZER0dTs2ZNAEaPHk1CQkKhX2vUqFFERPj/qzV16lSuvPJKGjZs\nSEREBDfddFNJlC0iImXcokVO5+cnn0CrVm5X446gCmvGmPrGmLnGmJXGmF+NMXdmHa9hjJljjFlt\njJltjAmjmerA++effxg8eDBdu3YlISGB7777DnDC2ty5cwv9enlt3v7uu++ydu1aevfuTbVwutBA\nRESKbOVKuOwymDIFunZ1uxr3BNs0aAZwr7V2uTGmCrDUGDMHGAR8a619xhjzIPAQMNzNQsuSP/74\nA4/Hw8CBA+nUqVNA3mPOnDnZ92fNmhWQ9xARkbIjJQX69IFx45x9P8NZUI2sWWs3W2uXZ93fD/wO\n1AcuBaZknTYFuMydCoPPn3/+Sd++falTpw4VK1YkLi6O/v374/F4ACeI9e3blxo1alCpUiU6derE\n7Nmzs58/aNAgzj77bAB69uxJZGQkgwYNIiIiAmMMjz/+OBEREURGRjJmzBhXPqOIiISXLVvg3HPh\n/vudKdBwF2wja9mMMQ2B1sAioI61dgs4gc4Yc4KLpQWVCy64gNjYWF577TViY2PZuHEjM2fOxOPx\nsGXLFrp06UK1atWYOHEiVatW5eWXX+bCCy/kq6++4rzzzmPEiBG0adOGu+66i1deeYUzzzyT2rVr\nc/vtt9OxY0cGDRrE4MGDAahfv77Ln1ZERMq63bvhvPOckHbnnW5XExyCMqxlTYF+DNxlrd1vjMm9\nr5D2GQJ27NjBmjVrGD9+PBdddFH28QEDBgAwbtw49uzZw+LFi2nUqBEA559/PqeeeiqPPPII5513\nHo0aNaJFixYAtGjRgvbt2wNkn1+vXr3sYyIiIoF08KCz32ePHjBihNvVBI+gC2vGmCicoPaOtXZG\n1uEtxpg61totxpgTga35PX/UqFHZ9+Pj44mPjy9OMUV/7vGUwL6WsbGxNG7cmOHDh7N582bi4+Np\n0qRJ9uOJiYl07NgxO3gBREREcPXVV/PYY4+xf/9+qlSpUuw6REREiuvIEbjiCmdpjvHjQ3Nj9uNJ\nSEgo0ioLQRfWgMnAb9baF3yOfQ7cCDwN3ADMyON5gH9YK7YQ2Cj822+/ZdSoUTz88MNs376dRo0a\n8cADDzB48GB27tzJWWedddRzTjzxRKy17Nq1S2FNRERcl5kJ118P0dHw5psQEVRX1Jec3INIo0eP\nLtDzgiqsGWO6ANcCvxpjfsKZ7nwYJ6R9aIy5CUgB+rlXZXBp2LAhb7/9NgC//PILL730EkOGDCEu\nLo6aNWuyefPmo57zzz//YIyhRo0apVytiIiIP2thyBDYtg1mzgyPjdkLK6iyq7V2obU20lrb2lp7\nprX2LGvt19bandbaXtbaZtba3tba3W7XGoxatWrFuHHjsNaycuVKevTowaJFi0hNTc0+x+Px8MEH\nH3DWWWcdd1StfPnyHDp0KNBli4hImLIWHnwQfvoJZswIn43ZC0v5NYT9+uuv3HXXXfTv358mTZqQ\nmZnJW2+9Rbly5TjnnHOoU6cOb7/9Nueeey6jRo0iJiaGiRMn8tdffzFz5ky/17J5TPmeeuqp2V2j\nNWrUoG7dupxUxA3Zfv/9d3777TestRw6dIiUlBQ++eQTwBkWjo2NLdLriohI6Bo7FmbNgoQEiIlx\nu5rgFVQja1I4J554InFxcYwfP55LL72Ua665hs2bN/PVV1/RunVrTjrpJBYuXMhpp53GkCFD6Nev\nH7t372bmzJmce+65fq+Ve8cBgJdffpnKlStzySWX0L59eyZNmlTg2nK/3ocffki/fv3o378/u3bt\nIiEhgX79+tGvXz9WrlxZtH8AIiISsiZMgLfecvb91N/Xj83kNaISqowx9lifxxiT5wiShB59lyIi\noWvKFHj0UZg/Hxo2dLsa92T9WXbcvldNg4qIiEip+eQTeOghmDs3vINaYSisSaFkZmYe8/HIyMhS\nqkRERELN7NlO5+fs2dC8udvVhA5dsyYFNm/ePMqVK5fvrXz58n6dpyIiIl7z5ztrqX32GbRu7XY1\noUXXrEmBHThwgNWrVx/znFatWhFVCovk6LsUEQkdixc720hNmwbnnON2NcGjoNesKaxJSNJ3KSIS\nGn75Bc4919mZwGcba6HgYU3ToCIiIhIQq1ZBnz7OMh0KakWnsCYiIiIlbt06Z0Rt7Fjop00ii0Vh\nTURERErUxo3QqxcMHw433uh2NaFPYU1ERERKzJYtThPB4MFwxx1uV1M2KKyJiIhIidi+3RlRu/pq\neOABt6spOxTWREREpNh274beveHCC2HECLerKVsU1qTENWzYkIEDB7pdhoiIlJJ9++D886FbN3jy\nSTDHXYxCCkNhTUqc0X+lIiJh4+BBZ1mOVq3g+ecV1AJBYU1ERESKJC0NLrsMGjWCV15RUAsUhbUy\n4Oeff6Zv377UqlWLSpUq0bx5c55++mkAvvnmGy688ELq1q1L5cqVadmyJc899xwej8fvNSIiIhgz\nZozfsZSUFCIiIpg6darf8Xnz5tG7d2+qV69OlSpVaN26NW+99dZRdX3wwQeceuqpVKlShXbt2rFw\n4cLsx5577jkqVKjAjh07jnpe48aNueaaa4r8z0NERALv8GG44gqoWdPZnSBCiSJg9I82xC1evJjO\nnTuzbt06XnjhBWbOnMl9993H33//DcDatWs5++yzeeONN5g5cyY33ngjo0eP5j//+U+R3m/GjBn0\n6tWLjIwMXn/9dT7//HNuvvlmUlJS/M6bP38+zz33HE888QQffvghmZmZXHzxxezduxeAQYMGERER\ncVTImz17NikpKdx+++1Fqk9ERALvyBFnodsKFeCddyAy0u2KyjbtDRriunfvTkpKCn/88QfR0dHH\nPT8zM5Onn36acePG+Y1qRUREMGrUKEb4tPCkpKTQqFEj3n777eyGgUaNGnHCCSfwww8/5PsejRo1\nYu/evaxbt46qVasCsHTpUtq1a8f777/PgAEDACewLVy4kD/++CP7uZdffjmrV69m5cqVx/wcZfG7\nFBEJBRkZMGCAE9g+/hjKl3e7otBV0L1Bo0qjmFBlEhIC9to2Pr7Yr3Ho0CGSkpJ48MEH8w1qmzdv\nZuTIkcyePZtNmzaRkZEBOP+CbN26lRNOOKHA77d69WpSUlJ4+OGHj3tup06dsoMaQMuWLQFITU3N\nPjZkyBCmTp3K3Llz6dmzJ5s3b+bLL7/kv//9b4FrEhGR0pOZCddfDwcOwPTpCmqlRWHtGEoiUAXS\nrl278Hg81KtXL8/HrbVcfPHFbN68mdGjR9OsWTMqVqzIZ599xtixY0lLSyvU+3lH4vJ7P181a9b0\n+7181n/Rvu/Zrl07zjrrLF599VV69uzJpEmTKFeunJb9EBEJQpmZMGiQs/Dt559DASZzpIQorIWw\nGjVqEBERwcaNG/N8fM2aNSxdupT33nuPq6++Ovv4jBkzjjo3OjqaI0eO+B3LffF/rVq1APJ9v6IY\nMmQIt912G5s2beLNN9+kX79+VK9evcReX0REis/jcbaP2rABvvoKKlZ0u6LwogaDEFaxYkW6du3K\nu+++y+HDh496/ODBgwBEReVk8vT0dN57772jzo2Li2PFihV+x7788ku/NdOaNm1Kw4YNeeONN0rq\nI3D11VdTpUoVrrnmGjZs2MDgwYNL7LVFRKT4PB647TZYvRq++AIqVXK7ovCjkbUQ99///pf4+Hg6\nduzIfffdR/369VmzZg2//PIL//3vf4mLi+ORRx4hIiKCqKgonn/+eSLy6K8eMGAATzzxBGPHjqVj\nx44kJiYybdq0o857/vnnueKKK+jZsye33XYbtWvX5vfff2fr1q2MGjWq0PVXqFCBG2+8kfHjx3PG\nGWfQsWPHovxjEBGRALDW2Yz9t99g1iyoUsXtisKTRtZCXNu2bVm4cCENGjTgzjvv5MILL2TcuHHU\nr1+fcuXKMWPGDE488URuuOEG/v3vf9OjRw+GDx9+1Os89NBDDB06lJdffpm+ffuyevVq3n333aPO\nu+SSS/jmm28wxnDLLbdw6aWXMmnSJBo1apR9jjEmz10M8jt+1VVXAWhUTUQkiFgL//43LF8OM2dC\nTIzbFYUvLd0hrnvkkUeYMGECmzZtokoB/9qm71JEJHCshbvvhkWLYM4cqFbN7YrKJi3dIUFv+fLl\nrFq1ihdffJHbbrutwEFNREQCx1q47z5ISoJvvlFQCwYaWRPXNGrUiK1bt9KnTx+mTp1K5cqVC/xc\nfZciIiXPWnjgAZg7F779FmrUcLuisq2gI2sKaxKS9F2KiJQsb1D77jsnqOVaLlMCoKBhTQ0GIiIi\nYU5BLbgprImIiIQxBbXgp7AmIiISphTUQoPCmoiISBhSUAsdWrpDREQkzHiX55g3T0EtFCisiYiI\nhBHvgrdJSVqeI1QorImIiIQJ7xZSS5Y4C95Wr+52RVIQumZNjislJYXRo0ezfv16t0sREZEi8nhg\nyBBYutTZQkpBLXQorMlxrV+/ntGjR7N27Vq3SxERkSLweOC22+CXX2D2bG0hFWoU1sLYkSNHCnSe\ntRZjjrvAsoiIBKHMTLj1Vvj9d/j6a6ha1e2KpLAU1kLcmjVrGDhwII0bN6ZSpUqccsopDBkyhN27\nd/udd+ONN3LyySezaNEiunTpQqVKlXjwwQezH580aRJt2rShUqVK1KxZk7PPPptFixYxb948evbs\nCUCvXr2IiIggMjKS+fPnA/DBBx9wzjnncMIJJxATE8NZZ53F1KlTj6ozIiKCESNGMGHCBBo3bkzV\nqlWJj4/nt99+C+A/HRGR8JaRATfeCGvXwqxZEBPjdkVSFNobNMQlJiYyc+ZMOnbsSM2aNVm3bh1j\nx44lNjaWhQsXZp83aNAgPvnkE2JjYxk2bBgtW7akYsWKtGvXjmHDhvHcc89x6623cskllxAREcGi\nRYs47bTTuOCCC3jnnXcYOnQoEyZMoG3btgCceuqpVKlShbFjx1K1alWaNm1KVFQU8+fP58knn2TC\nhAn861//yn7/iIgIGjZsSLNmzbjjjjs4cuQIw4YNIyoqilWrVhERUbi/N5TF71JEpCSlp8PAgbB9\nO8yYAZUquV2R5FbQvUGx1paZm/Nx8ne8x8uCjIwMu2DBAhsREWGXL1+effzGG2+0ERER9osvvvA7\n/6+//rKRkZF22LBh+b5mQkKCNcbY77777pjv7fF4bEZGhr311ltt69at/R4zxtimTZvajIyM7GMf\nf/yxjYiIsMnJyYX5iNba8PguRUSK6vBha6+4wtrzz7f20CG3q5H8ZP1Zdtx8o6U7jiHBJATsteNt\nfIm8Tnp6Os8++yzvvPMOKSkppKWlAU5aX716NWeccUb2ueXKlePCCy/0e/63336LtZZbb721SO//\n119/8eijj5KYmMjmzZvxeDwAVKhQ4ahzzz33XCIjI7N/b9myJdZaUlNT6dixY5HeX0RE/B0+DP36\nOfc/+wyio92tR4pPYe0YSipQBdLw4cN5+eWXGTlyJJ06dSImJoa///6bvn37Zgc3r9q1ax/VKLBj\nxw4A6tevX+j3PnDgAL169aJKlSo888wzNG7cmPLlyzNx4kTeeuuto86vmWuJ7Ois/4PkrlNERIom\nLQ0uv9yZ8nz/fShf3u2KpCQorIW4Dz74gBtuuIGHHnoo+9i+ffvyPDevjs5atWoBsHHjRv7v//4v\n3/fJ67nJycls2LCBBQsW0KlTp+zj6enpBa5fRERKxoEDcOmlcMIJMHUqROlP+DJD3aAh7uDBg0Tl\n+i9y8uTJBV5qo1evXhhjeP311/M9Jzo6Gmsthw4dOuq9Ab/337VrF59//nlByxcRkRKwdy/06QMN\nGsA77yiolTX6OkNcnz59mDJlCqeffjpNmjTh008/JTk5ucDPb9y4Mffccw/jx49n7969XHLJJURG\nRrJ48WJatGjBVVddld3pOXnyZGrUqEF0dDTNmzenc+fOxMTEcMcddzBq1Cj279/PE088Qe3atdm7\nd28AP7WIiHjt3OkEtbZt4aWXoJDN9RIC9JWGuAkTJnDJJZfwn//8hwEDBnDgwAH+97//5XlufqNt\nzz77LBMnTuSHH37gyiuv5LrrriMhIYEGDRoAzrVmL7/8Mj///DPx8fG0b9+epUuXUqtWLaZPn05m\nZiZXXXUVjzzyCLfeeivXXnttnu+d1/trsV0RkaLbtg169oRu3eDllxXUyiqtsyYhSd+liIS7f/6B\nXr2choIxY0B/9w09BV1nTRlcREQkxKSmQo8ecM018NhjCmplncKaiIhICPnzT+jeHYYMgUcecbsa\nKQ1qMBAREQkRK1bAeefB6NFwyy1uVyOlRWFNREQkBCxdChdeCOPHw9VXu12NlCaFNRERkSC3YIHT\nSDBpkrPwrYQXhTUREZEg9s03TiPB++/Duee6XY24QQ0GIiIiQeqzz+Daa+HTTxXUwpnCmoiISBCa\nOtXp+PzUvCWKAAAgAElEQVT6a2fRWwlfQTcNaox5E7gI2GKtbZV1bCRwK7A167SHrbVfF/a14+Li\ntGJ+GREXF+d2CSIiAfPSS/D00zB3LrRo4XY14rag28HAGNMV2A9MzRXW9llrnzvOc4+5g4GIiEgw\nsxbGjoW33oJvv4WGDd2uSAKpoDsYBN3ImrV2gTEmr2ETDYmJiEiZZS088IAz7ZmYCCed5HZFEixC\n6Zq1ocaY5caYN4wx1dwuRkREpKRkZsK//gXz58O8eQpq4i9UwtpEoLG1tjWwGTjmdKiIiEioOHwY\n+veHdevgu++gZk23K5JgE3TToHmx1m7z+XUS8EV+544aNSr7fnx8PPHx8QGrS0REpDj274e+faFa\nNfjqK4iOdrsiCaSEhAQSEhIK/bygazAAMMY0BL6w1rbM+v1Ea+3mrPv3AO2stdfk8Tw1GIiISEjY\nscPZPur00+G11yAy0u2KpLSFbIOBMeZ9IB6INcakAiOBs40xrQEPsB4Y7FqBIiIixbRxI/TuDRdd\nBE89BVpVSo4lKEfWikojayIiEuz+/BPOOw8GD4YHH3S7GnFTyI6siYiIlFXLljmjaY89Bjff7HY1\nEioU1kRERErB3LkwYAC8/jpcdpnb1UgoCZWlO0RERELWJ584Qe2jjxTUpPAU1kRERALotdfg3/+G\n2bOhRw+3q5FQpGlQERGRALAWHn8c3n7b2ZmgSRO3K5JQpbAmIiJSwjIz4c47YeFCWLBA20dJ8Sis\niYiIlKC0NLjuOti509nns5p2s5Zi0jVrIiIiJWTPHujTByIiYNYsBTUpGQprIiIiJWDTJujeHVq2\nhGnTtM+nlByFNRERkWJavRq6dIH+/eHFF7XPp5QsXbMmIiJSDMnJ0LcvjB0LN93kdjVSFimsiYiI\nFNGMGXDLLTB1Kpx/vtvVSFmlaVAREZEieO01uP12mDlTQU0CSyNrIiIihWAtjBgB//sfJCbCKae4\nXZGUdUUOa8aYhkAPoDlQAzgIbAV+BuZaaw+XQH0iIiJBIz0dBg+GFSucBW9POMHtiiQcFHoa1Bhz\ngTEmEfgfcCawB1gKrAWigauBX4wxrxlj6pRksSIiIm7Zuxcuugi2bYPvv1dQk9JT4JE1Y0x5YByw\nE7jUWrvzOOe3B543xnxlrX23eGWKiIi4Z9MmuOAC6NQJJkyAKF1EJKXIWGsLdqIxTwJvWGvXFOoN\njBkI7LbWfl6E+grFGGML+nlEREQKYuVKJ6jdfjs8+CAY43ZFUlYYY7DWHvffqMKEtfLW2iNFLKbI\nzy3k+yisiYhIifn+exgwAJ57Dq691u1qpKwp8bAWChTWRESkpLz/Ptx9t9P12bOn29VIWVTQsKZZ\ndxERER/WwhNPwKRJMHcunH662xVJuDtuN6gx5gxjzM/GmN3GmP8ZY2pnHb/GGDMr8CWKiIiUjvR0\nuPlm+OwzWLRIQU2Cw3GnQY0xnwFTgD+BzsBVwEBr7WZjzFZrbdA0L2saVEREimrPHrjySqhQAaZN\ngypV3K5IyrqCToMWZJ21L6210621K621k4Argbu1hpqIiJQVKSnQpQs0bw7TpyuoSXApSFizxpjT\njTETjDHVrLV7gYeAS4AKgS1PREQksH780QlqN98ML74IkZFuVyTi77gNBtbaycaYPjjToPuzjllg\nkjFme4DrExERCZjp0+HWW51mgssuc7sakbxp6Q4REQk71sK4cfD88zBjBrRp43ZFEo5KZekOY8wY\n4HzgSmttSnFeS0REpDSkp8PQoU63Z3IynHyy2xWJHFtx11nbDcwEPCVQi4iISEDt2QNXXeXs7blg\nAcTEuF2RyPEVpMHgWHYBT1trN5REMSIiIoGydi107gxNm8LnnyuoSegoblhbAnxpjLneGFO/JAoS\nEREpaQsWOB2ft98OL73kjKyJhIpiNRgYY94FNgAdgU5AKjAX+NRaO6dEKixcPWowEBERP++8A/fd\nB1OnQp8+blcjkqNUNnI3xjwEzLDW/maMqQh0BXoCzay1lxf5hYtej8KaiIgA4PHAo486uxF88QWc\ndprbFYn4K5WwlvVGvYGa1tr/FeuFSoDCmoiIABw8CAMHwubNzj6ftWu7XZHI0UpyuynvC9bK67i1\nds7xgpp383cREZFA+/tv6NYNKlWC775TUJPQV5gGg0bGmFsK+wbGmI7A7YV9noiISGH98AN06AD9\n+sGUKRAd7XZFIsVX4H4Ya+0SY4w1xnwK/A+niSAjv/ONMa2AO4EUa+2Y4pcqIiKSv/ffh7vvhjfe\ngEsucbsakZJT6GvWjDHRwF3A9cAOYDXO4rhHgJrAiUArYDEwylr7Z0kWfJzadM2aiEiY8XhgxAh4\n7z1n/bSWLd2uSKRgSqsbtBlwJlAHiAa2AuuAJGttepFfuOj1KKyJiISR/fudRoJt2+DTT3V9moSW\nUusGzXqzWGCn20lJYU1EJHysWweXXgrt2sHEibo+TUJPiXeD5vMm8caYDTgjajuMMc8ZY6oU5zVF\nRESOJyEBOnWCW291rlFTUJOyrLgbbgwELsq63w7oBywyxnSz1u4q5muLiIgc5ZVXYNQop6HgnHPc\nrkYk8Ip7zdqd1toXcx3rA5xvrb2ruMUVoR5Ng4qIlFFHjsCdd0JiotNIcMopblckUjylMg0K1DPG\n+G3gYa39Gme/UBERkRKxdSv06gWbNkFysoKahJfihrVngClZ16q18Tle3OlVERERAJYudZoIzj4b\npk+HqlXdrkikdBUrrFlrdwBnAx7gW2PMdmPM70ANY0zjkihQRETC1zvvQJ8+MH48jB4NEcUdYhAJ\nQSWydAeAMSYS6Aj0BM7Jur8R+MBa+3CJvMnxa9A1ayIiZUBGBjzwgHNt2vTpcPrpblckUvJKdZ21\nfAqoCHQHTrfWjgvImxz9ngprIiIhbvt2GDAAIiPhf/+DGjXcrkgkMEqrwSBf1tpD1trZpRXUREQk\n9C1b5lyf1qYNzJypoCYCagQQEZEg8c47cO+98PLL0K+f29WIBA+FNRERcVV6Otx3H8yaBd9/r+vT\nRHJTWBMREdds2QJXXQUxMbBkCVSv7nZFIsFHTdAiIuKK5GRo2xbi4+GLLxTURPKjkTURESlV1ubs\n7/nmm3DxxW5XJBLcFNZERKTUHDwIt90Gy5dDUhI0aeJ2RSLBT9OgIiJSKtasgU6dwOOBRYsU1EQK\nSmFNREQC7ssvoXNn+Ne/nCU6KlVyuyKR0KFpUBERCZjMTBgxAqZOhc8+cwKbiBRO0IU1Y8ybwEXA\nFmttq6xjNYAPgDhgPdDPWrvHtSJFROS4tmyBa64BY2DpUjjhBLcrEglNwTgN+hZwXq5jw4FvrbXN\ngLnAQ6VelYiIFNjChc6yHJ07w+zZCmoixRGwjdyLwxgTB3zhM7K2Cuhhrd1ijDkRSLDWNs/jedrI\nXUTERdbC88/DU0/BW2/BBRe4XZFI8CroRu5BNw2ajxOstVsArLWbjTH6O5qISJDZvRtuugk2bIAf\nfoCGDd2uSKRsCMZp0ILQ8JmISBBZtgzatIF69WDBAgU1kZIUKiNrW4wxdXymQbfmd+KoUaOy78fH\nxxMfHx/46kREwpS18Npr8Oij8PLL0K+f2xWJBK+EhAQSEhIK/bxgvWatIc41ay2zfn8a2GmtfdoY\n8yBQw1o7PI/n6Zo1EZFSsm8fDB4MK1fCRx9B06ZuVyQSWgp6zVrQTYMaY94HkoCmxphUY8wg4Cng\nXGPMauCcrN9FRMQlP//sdHtWquTsRqCgJhI4QTmyVlQaWRMRCSxr4fXX4T//gfHj4brr3K5IJHSV\ntW5QERFx2d69znZRv//uNBE0a+Z2RSLhIeimQUVEJPgsWwZnnQXVqzvTngpqIqVHYU1ERPJlLUyY\nAOedB088Aa++ChUrul2VSHjRNKiIiORpxw5nkduNGyE5GZo0cbsikfCkkTURETlKYiKceSb83/9B\nUpKCmoibNLImIiLZMjNh7FiYOBHefFN7e4oEA4U1EREB4O+/4frrwRhYuhTq1nW7IhEBTYOKiAgw\nfbqzt2evXvDNNwpqIsFEI2siImHs4EG4916YMwdmzICOHd2uSCQ87Du8r8DnKqyJiISpX36Bq6+G\n1q3hp5+gWjW3KxIpu7Yd2MaC1AUkpiYyP2U+q7avKvBztd2UiEiY8XictdMefxyee865Tk1ESlbq\nnlQSU5xglpiayKZ9m+h8cme6NehGt7hutKvbjgrlKhRouymFNRGRMPLPPzBoEOzeDe++qyU5REqC\ntZbVO1aTmJKYPXJ2MP0g3eO6Z4ezM+qcQWREpN/zCro3qMKaiEiYmDEDBg929vd89FEoV87tikRC\nU6Ynk5+3/Jw9apaYkkilcpXoFteNbg260SOuB01jm2LMsXOYwpqIiABw4ADcdx/Mnu2MpnXp4nZF\nIqElLSONJRuXOMEsNZGkDUnUi6nnN3LWoFqDQr+uwpqIiPDjj3DttdChg3OdmpoIRI5v3+F9JG1I\nyp7SXPbPMlrUbkG3Bt3oHtedrg26UqtSrWK/j8KaiEgYy8iAJ590AtoLLzhdnyKSt7w6NdvUbeOM\nmjXoRueTOxMTHVPi76uwJiISpv76y+nwrFwZ3n4b6td3uyKR4OLt1PSGs437NmZ3anaP607bum2p\nEFUh4HUorImIhBlr4Y034KGHnAaCf/8bIrRPjYQ5b6embzPAwfSD2c0A3eO606pOK6IiSn/pWYU1\nEZEwsnkz3Hqrs7/ne+/Bqae6XZGIO/Lr1PRtBmgW2+y4nZqloaBhTTsYiIiEuI8/hqFD4eab4ZNP\noHx5tysSKT2HMw6zZNMSZwHa1PkkbUiiftX6dGvQjcubX87488YXqVMzoHbuhKSkAp+ukTURkRC1\na5cz1bl4MUydqn09JTz4dmompiaydNPSgHRqlhhrITUVEhNhwQLnlpoKHTpgvv1W06AiImXV7Nlw\nyy1w2WXw1FNOM4FIWXSsTs3ucd3pVL9TQDo1iywzE1asyAlmCxZAejp07Qrdujk/zzgDoqJ0zZqI\nSFm0fz/cfz989RVMngy9erldkUjJSt2T6lxvltWt6VanZoEdOgRLluQEs6QkqFPHP5ydcgrkcY2c\nwpqISBmTkAA33QTx8c4G7NWru12RSPHk7tScnzKfQ+mHjrunpqt27oSFC3PC2fLlcNppTijr2tXZ\nIqROnQK9lMKaiEgZceCAsxzHp5/Ca6/BhRe6XZFI0RxrT83uDbrTPa57gfbULDXWQkqK/5Rmaqpz\ngWiXLs7IWYcORb4OQWFNRKQMWLAAbrwROnWCF1+EGjXcrkik4Lydmt5wlrwhmboxdYu9p2bAeK83\n820GyMjImc70ud6sJCisiYiEsAMH4JFH4MMPYeJEp5FAJNgdq1OzW4NudG3QldqVa7tdZo5Dh5x2\nam8wS06GE0/0D2eNG+d5vVlJUFgTEQlR8+Y5a6Z17Ojs6xkb63ZFInnzdmp6R85Ka0/NItu+3f96\ns19+gZYt/a83q116YVJhTUQkxOzfD8OHw/TpzmjaJZe4XZGIvw17NjA/ZX52OPN2anZv0J1ucd1o\nV7cd0VHRbpfpsBbWrXNCmXdac9Omo683q1TJtRIV1kREQsh33znrpnk7PXVtmrjN26np3RnAd09N\nbzgLqk7NjAxnpMy3GcAY/yU0WraEyCCpF4U1EZGQsHu3s27a7Nnw6qtwwQVuVyTh6nidmsG0pybg\nXNj5ww85wWzRIjj55Jwpza5doWHDgF1vVhIU1kREgtz06XDHHXDppc4uBFWrul2RhJO0jDSWbFyS\n3QyQtCGJejH1grdTc8sW/+vNVq6E1q1zglnnziF3gafCmohIkNqyxdnTc/lyeOMN6N7d7YokHOw7\nvI/kv5Ozrzlb9s8ymtdqTve47sG3p6a18Oef/lOa27Y5gcx7vVnbtlCxotuVFovCmohIkLHW2XD9\ngQdg0CAYOTLk/6yRIOa7p2ZiaiK/b/udNnXbZE9pBtWemunp8NNP/uGsYsWcYNa1q7NLQESE25WW\nKIU1EZEgsmYN3HYb7NgBkyZBmzZuVyRlTeqe1Oz9NOenzA/uTs29e51rzLzBbMkSZ/9M7/IZXbpA\ngyCagg0QhTURkSCQng7jx8Mzz8CDD8I995TY4ucSxnLvqenbqend8DyoOjU3bnRCmfeasz/+cKYx\nu3RxAlqnTmG52a3CmoiIy3780VmOo3ZtZ0/Pxo3drkhCVYYng583/5w9pent1PRtBgiaTk2PB377\nzb8ZYN++nFGzrl2doeXy5d2u1HUKayIiLtm3Dx59FKZNg2efheuvD+rVAyQI+XZqzk+ZT/LfydSv\nWj971Kxbg26cXO1kt8t0pKU505jecJaU5HRleoNZ167QrJn+I8iDwpqIiAumT3c6PXv1coJarSBp\nrpPgdqw9NbvHdafLyV2CZ0/NHTucQOYdNVu+HE491X/LphNPdLvKkKCwJiJSilJTnZC2erWzuG18\nvNsVSTDLb0/NoOvUtBbWr/fv0tywwdmmyRvOOnSAKlXcrjQkKayJiJSCjAx48UUYOxbuvNNpIogO\nkoY7CR7eTs2g31Mz95ZNCxc6gc07YtatG7RqpS6ZEqKwJiISYElJcPvtTgPBxInQtKnbFUkwOFan\nZtDtqbl/v7OEhvd6sx9+cLZs8l3fLMi3bAplCmsiIgGyYwcMHw5ffQXjxsGAAfqzLJyF1J6amzb5\nd2muWgVnnum/ZVPNmm5XGTYU1kRESpjHA2+/DQ89BP37w2OPQbVqblclpe1wxmGWbFqSHc68e2pm\nd2oGy56aHo8TxnynNHfvzll01ruERoUKblcathTWRERK0PLlzqbr6enwyivagSCc+HZq5t5Ts1uD\nbnRt0DU4OjXT0pzF/bzBLCnJWWjWt0uzefMyt2VTKFNYExEpAbt3O2umffghPP443Hyz/qwr6/Lr\n1PSOnAVNp2buJTR+/hlatMgZNevSBU46ye0q5RgU1kREisHjcTZdf+ghuPRSeOIJZ51PKXtColPT\nWmeDWe/1ZgsXOls4deyYE87at9cSGiFGYU1EpIiWLXPWTEtPd7o827Z1uyIpKd5OzcSUROanzs/u\n1PTdtikoOjXT0525d9/9NKOi/Kc0W7bUEhohTmFNRKSQtm+HRx6BGTOc5gFNeYa+/Do1veGse1x3\nmsY2db9Tc88eSE7OCWY//giNGuUEsy5dIC5ObcdljMKaiEgBZWQ4G62PHg1XXw2jRkGNGm5XJUVx\nrD01vSNnQdGpuWGDf5fmX385Q7jeKc1OnZzmACnTFNZERApg3jxn54HYWHjhBWdmSUKHt1PTO3K2\n7J9l2XtqBk2nZmYm/Pqr//VmaWk5o2ZduzprnZUv726dUuoU1kREjmH9erj/fliyxNlw/corNcMU\nCrydmt6Rs6DcU3P/fmcngIULnduiRVC3rn+XZpMm+hdOFNZERPKyfz889ZSzVto998B990HFim5X\nJfnxdmp6w1lQdmr67gqwcCH8/ju0bp0TzDp3hlq13K1RgpLCmoiID48H3nvPWYrj7LPhySehfn23\nqxJfvp2a3nDm3VPT2wzQqk4roiJc7ID0eOC33/zD2Z49TiDzjpy1batdAaRAFNZERLIsWAD33uvc\nf+EF59ptcZ+3U9O7jMaC1AVUjKrot4yG63tqHjzozJV7w1lysjNK5jul2ayZ2oalSBTWRCTsrV0L\nDz7oXD701FPOhuv6M9U93j01veEseUMydWPq+oUz1zs1N2/OudZs4UJYscLpOvGGs86doU4dd2uU\nMkNhTUTC1p49zo4Dkyc716Xdcw9UquR2VeHHd0/NxNRElm5amr2nZve47nRt0JValVy8lsvjca4v\n8w1nO3bkTGl26QLt2ulfHgkYhTURCTtHjjjrpT3xBFx8sbOw7Yknul1V+MivUzNo9tQ8dChnStO7\n0XmNGjnBrGtXZ29NDb9KKSmTYc0Ysx7YA3iAdGtt+1yPK6yJhCFr4dNPYfhwOOUUeOYZaNXK7arK\nvvz21PSGM9c7Nbds8R81+/VXOP30nHDWubM2OhdXldWwthZoY63dlc/jCmsiYSY5GYYNgwMHnPXS\nzj3X7YrKJm+npjeYzU+Zz6H0Q3SL65a9jIare2p6PLBqVU4wW7DAmdLs1Cln1ExTmhJkympYWwe0\ntdbuyOdxhTWRMPHHH/Dww07zwOOPw3XXQaTLe2+XJUG/p2buKc3kZGd7Ju+oWZcucOqpmtKUoFZW\nw9paYDeQCbxurZ2U63GFNZEy7p9/YMwY+OgjZ0Ttzjs1WFIS8tpTs15MveDp1MyrS1NTmhLiympY\nO8la+48xpjbwDTDUWrvA53GFNZEyau9eZ5pz4kQYNMhZ3DY21u2qQpdvp+b8lPnBtaem78Kz3tuu\nXTlTmurSlDKioGHNxWWgC89a+0/Wz23GmM+A9sAC33NGjRqVfT8+Pp74+PhSrFBESlpamrM11FNP\nwfnnw7JlEBfndlWhx9up6Z3WXLV9FWeddBbd47rzn+7/cbdT88ABWLzYfy9N78Kz3bo5nSPNm2tK\nU0JeQkICCQkJhX5eyIysGWMqARHW2v3GmMrAHGC0tXaOzzkaWRMpIzIyYMoUGD0azjzTuS6tZUu3\nqwodqXtSnWCWtXVTUHVqbtzoP2r2++9wxhk565tp4VkJE2VuGtQY0wj4DLA4I4LvWWufynWOwppI\niPN44OOP4dFHoW5dGDtW20MdT0E6NV3bUzMzE375xVnTzBvODhzwD2bt2mkvTQlLZS6sFYTCmkjo\nshZmzXJCmjFOSDv3XOe++MuvU9M3nLm2p+bevc40pjec/fAD1KuXE8y6dIGmTfXFiqCwJiIhZO5c\n+M9/nD/nR4+Gyy/Xn+W+grZT01pYv95/1GzNGjjrrJxGgE6d1Akikg+FNREJegsXOiHt77+dkNa/\nv9ZKg5xOTe/ImW+npqt7ah45Aj/9lBPOkpKcwOa7tlnr1lC+fOnXJhKCFNZEJGgtWgSjRjkLzo8Y\nAQMHQlRI9aaXrLw6NdvUbZM9pelap+aOHU4g84azZcugSRNnOtN7a9RIw6AiRaSwJiJBZ/FiGDkS\nVq6ERx5x1ksLx0GY/PbU9IYzVzo1PR5Yvdo/nP3zD3TokDOd2bEjVK1aunWJlGEKayISNJYscUbS\nfv3V2SJq0CCIdnF/79Lk7dRMTElkfqqzlMbB9IPu76l54IDzxXjDmXe7Jt+FZ087TfPSIgGksCYi\nrktOhscec0LaQw/BzTeX/ZB2vD01XevU3LDBf9Ts99+dhet8GwG0XZNIqVJYExHXJCY6Ie2PP5zF\n58vySFpenZr1q9bP3rbJlU7N9HT4+Wf/RoDDh3OuM+vSBdq00dpmIi5TWBORUmUtfP+9E9JSU53p\nzuuvL3vXpOW1p2bzWs3pHtfdvU7NHTucYUzvyNnSpc6F/95w1qmT0xigRgCRoKKwJiKlwlr48ktn\nEdudO52Qdu21Zae709up6Q1n3k5N7zIapd6p6fE4bbTeYJaUlNMI4L3erEMHqFat9GoSkSJRWBOR\ngMrMhI8+giefdAZsHnnEWcw21K9Hz69T0xvO2tZtS4WoUpw+3L/faaP1bQSIjc0JZp07qxFAJEQp\nrIlIQBw+DO+8A888A7VrOyHt/PNDc4btWHtqesNZqe6paS2sW+c/pfnnn85Cs75TmtrkXKRMUFgT\nkRK1dy+89ho8/zy0agUPPgg9eoRWSDvenprd47rTNLZp6XVqpqU515d5R8ySkiAiImfErHNnOPPM\nsnfhn4gACmsiUkK2bIEXXoDXX4feveGBB5yBnlDg26mZmJpI0oak7E5N71IaJ1c7ufQK2rjRP5j9\n+iu0aJEzYta5MzRoEFoJWESKTGFNRIpl9WoYNw4+/hiuvhruuw8aN3a7qmMLqj0109Nh+XL/cHbw\noBPKvNebtW0LlSuXTj0iEnQU1kSk0KyFBQvgv/919u8cMsS51a7tdmV5y72n5u/bf6dt3bbZa5x1\nPrlz6XVqbtnihDJvMPvpJyfdekfMOnfW8hki4kdhTUQKLCMDPvvMGUnbscMZRRs4ECpVcrsyf0HT\nqZmRAb/8khPOkpOddUs6dswZOevQQftoisgxKayJyHHt2QNvvgkvvgj168O998KllwbHKhC+e2p6\nOzUPph/MvtbM26lZKntqbtvmH8yWLnWuLfMGs06doHlzpzlARKSAFNZEJF9r1zoBbepUOO88uOce\naN/e3Zpyd2ouSF1AxaiKfuGsVDo1MzJgxQr/cLZ1a86is506OSNo1asHtg4RKfMU1kTEj7WQkAAT\nJsD8+c6m6kOHwsml2Azp63DGYZZsWpIdzpI2JFEvpl5Op2Zp7am5fbtzgZ43mP34I9Sr5z9q1qJF\ncAw3ikjISfd4OJCZyf48bhfXrq2wJiJOA+L77zsjaRkZcOedcN11UKVK6dYRFHtqZmQ4y2UkJ+cE\ntG3bnGFFbzBr3x5q1gxsHSISlNI9Hr8wtS9XuDrq94yMfB/3BrR0a6kcGUlMZCRVct1mnnGGwppI\nOFu/Hl55BSZPdmbt7roLzjmn9JoRc3dqurKn5tat/sFs6VJnKNF3OlOjZiIhyVpLWla42ue9ZWRk\nB6Z9PuEp9/H9uR7zBqwMa6kSGUlMVJTzMzIyO2jlFbbyOpY7mEVHROR7+YamQUXCkMcDc+bAxInO\n6hEDBzpLbzRpEvj39nZqekfOvJ2a3Rs4U5rt6rYjOio6cAWkp8PPP+cEs0WLnA7NDh1yujTbt4ca\nNQJXg4gcU4bHkx2i9mYFqNxBK8/H8whh+zIyiDKGmKio7HDkG6p8A5ff4z7HcweuCscIVoGgsCYS\nRnbtgrffdkJalSpwxx3OQraBWm/Vt1Nzfup8ElMSOZh+MHvbpm5x3TijzhmB7dTcuNEJZN5wtny5\ns65Zx445zQDq0BQptkxr2ZeRwV6fELU3Kyzt9QlVvsf25brvPeewx5MdmKr6hKcYn1BVmN/Lhfh/\n369AOXgAABflSURBVAprImWctbB4Mbz6qrNG2gUXOCGtc+eSn+o83p6a3eK60Sy2WeD+RpqWBsuW\n+YeztDT/YNaundY1E/FxxOPJDla+P/fk+v14P9NyBayqWWHJ777PsdxBzPd+pVIeuQp2CmsiZdT+\n/U7DwKuvOuukDR4MN94IJ5xQcu/hu6fm/JT5JP+dTL2YetnLaAS0U9NaWLcuJ5gtWgQrVzqjZN7r\nzDp2hFNO0W4AUib5jmTtychgj8/93Mf25r7vE7QyrKWaT5Cq5hO2qhbieOXISAWsAFFYEylDrHWu\njZ80CT76COLj4bbboFevkpnlO1anZrcG3ejaoCu1Kwdoz6m9e2HJkpxg9sMPUK6cfzA766zg205B\nJA/eoLU7Kzh5g9XurJ++x3x/3+tz/0BmJpV9wpLvT+/NG6qqeR/Pul/V59zSvv5KCk9hTaQM2LMH\n3nvPCWl79sAttzijaHXrFu91vZ2a3nC2avsqzjrprJxOzZM7UTU6AFOKmZnOKNkPP+QEs/Xr4cwz\nc4JZhw7OdgoiLkjLzMwOWruzQtbuXIHrqGNZIWt3VtCq4hOsqvsEquzbMX73TiNGKmSFBYU1kRBl\nrbNo7eTJ8Pnn0Ls33Hor9OxZ9FE01/bU3LTJCWTe29KlcNJJOaGsQwdo1coZSRMpAUc8nuwwtcsn\nWO1KT/cLX/ndPED1rJBV3SdwVc8VvHwf9w1lClpSGAprIiFmwwaYMsXp6qxQAW66Ca6/HmoXcvax\nIJ2areq0IioiqmQ/wMGDThjzDWcHDuSEso4dtXSGHJe1lv2ZmdlBa5dP0NrlE8ByH/eGrSPWUj0q\nihpZAapG7rCVxzHf45o6lNKksCYSAg4dghkznIC2eDH07++EtLZtC37tvLdT0zecBbxT0+OBVav8\ng9nq1XDaaf6jZk2aqAkgDFlrOeTxsCsjg53p6f4/s0LWLp/wtdPn990ZGZQ3hhpRUdQoVy47cNXw\nDWDlyh0Vxrw/dTG8hBKFNZEgZS0sXOiMon3yiRPMbrgB+vYt2DX0vp2a3j0161etT9eTuwZuT83N\nm3NC2eLFzv6ZsbE5oaxDB2jd2hkSlDLDYy17swLWzvT07J+78jjmDWHe3wFqZoUq78/c92v4BC/v\nY9WjoogO8bWzRApKYU0kyKxZA+++C1OnQnS0E9Cuvfb419Ln1anZonaL7OvNSnxPzf37nenMxYud\nm3c6s317J5S1b+/cagV4H08pMdZa9mVmsiMrTO1IT/e7v9P3fq7wVTkyMjtgxZYrR82sUFUzK2DF\n+gSymj5hrKK28BI5LoU1kSCwbRt8+KET0tascaY5Bw489jRnqe6pmZEBK1bkBLPFi51CW7XKCWft\n2mk6M4ikezx+gWtHejrb09PZkeuY7++7MjKIjogg1hu4soJVbLly2QHMN4h571ePigr5FeJFgpnC\nmohLDhxwujjfew8WLIALL3RG0M49N++mx2N1anZr0I129dqVTKemd7FZ32C2fDk0aJAzWta+vRPU\nypcv/vvJcaV7PNlhyxu4tvv+nsfP/ZmZ2YEq1idweX+v5fO7byArr9AlEnQU1kRK0eHD8PXXMG2a\n87NTJyegXXaZs1enl7dT03fbpoDtqblli7PY7JIlTjBbssS5pqx9e2e0rEMHaNMGqlUr/nsJ1lr2\nZGSwLStYbfMJXb73fW/7MjKoWa5cdsCqnXXf+3stnwDmPVYtKur/27uz4DiOw4zj/957FzeJGyAh\nURJ1UZRIiYqsg5HkSFHkiu0qHyVXqpzyQ56ScioPqRwPSR6TpxxPebCTipOKHVu2RcuxLFuWVIRM\nheIFXuJlSiIJAiAAgiAJYHdndqbzMLPAAAQpkBK5u8D3q5qa7p7ZYVNLCh97enqIaZRTZFlQWBO5\nyVwXfvWr4DbnK6/Agw/CSy/Bl740N53rau/ULL8ZYGvfVtavXv/Jn167dGlunlk5oF26FISycjjb\nsuWTr6a7gvjWMhEGrTHXZdRxZstjrsvYgvq465KNxWgNQ1c5eLWlUrNhKxrGWsPJ9ApeIiuXwprI\nTVAOaD/4QbDkxvr18JWvwFe/Cj09t+idmvl8cPuyHMp27YLBwSAtRsOZ5pnNUx75Gg2D12gkgM0r\nh/uJUomGeJz2MHC1JZO0RwJY24JQ1ppMktGkehG5DgprIp+SYhHeeCNYZmPbNrj77iCgffnL0Nx+\nk9+p6bpw8GCwVMauXcH+2LHgpebl0bItW4L1zRKf8iK3NcDx/XnB61wkiJ1znCtCWToWoyOVoj2Z\npD3ct4XlefvwlqMm14vIzaSwJvIJTE0Fc89+9CN47TV44IFgHbRnPjfGh96V79Qsh7NP9E7NUgmO\nHAkC2Z49QTg7dAhuvz14fHTLlmD/4IPLej2zgudxLgxb5W0kDGLRtlHX5bLn0ZZMzgtgHeHIV0e5\nPRLKNPIlItVEYU3kOo2NwauvBqNnb70VPCTwzBdP07ihnwMXgwVoBy8N8viax2cfBtjSvYV0In39\nv5jnBSNk5WC2ezfs3x/cS33kkWB7+GHYvHn+Ewo1yvV9xlyXkTB4LdzORfYzvk9HJGx1plJBPdJW\n3lo050tEapjCmsgSnDgRhLNt2+DAQctjnztG7xP9TK3azs6RT+lJzXIw27NnbhsYgI6OuVBW3mro\nyUxrLZc8j+FikeEwbA1fJYxdKJVoTSbpDMNXeesI26KhrCWR0OuCRGRFUFgTWYTrwo4d8NOfwqv/\n6zEe38/657ZjbuvnRHH+OzVv6ElNzwvembl375XBLBrKNm+u2hea+9Zy3nUZDsPXUBjGhh1nNpiV\nQ1ncGLpSKbrCsNWVStGVTl8RylqTSeIKYCLLhrUWLGAjZZhtW3hO+djsz2h7lXOv9xpc49jC61zl\nM1ccX3i9pXzmWv2I7heUGx9pVFgTgeD25uuvw09+VuC1A7to3NBP7p5+hpM7WNv8CZ7UdF14//25\nYLZ3Lxw4ECyPsWlTEMoeeSQoV0Ews9YyUSoxVCxyNgxdQ5EwNlQsMhTeiqyPx2eDVzmMdaVSdEfq\nnakU9SvwoYZbyVqL9Sx4YD2L9cOyH2lfrOwzb2/9BW3R8yLXXPTYje4ti7bd8Dk2aJs9bpd4vNy2\nyPEr2q9WXnh+9NeMhIJ5AcZf5Jh/5Xm3tM4NnBOpzwscZWZum/2HbbmNSLvhimOz7QvPvd5rwJKv\nf9XPLDy+4HpL+sy1+hE9Filv2btFYU1WJs8L5ua/8tplfrxrBx/Zfho2bGeqYS/3rL6XZ++8gXdq\n5vPBU5n79gWhbO9eOHwYbrstGCUrj5Y99FBFbmVOex5nwxBWDmDzyuGoWC4epzsMYd1h+OqOhLDu\nMIRVy0R8ay3WDTbf8Rcvu3Nl61psKdJWWqS9NFe/onytzbt6HY/5bZH9bNAqt0fbPqaOBWJg4gYT\nNx9bXvR4bJHjC9tic5/BLLhG+VxjIB6px8zsufPOK7fH5rfP/sCLLzgnco15nzPzrzfbVi4vuEb0\nM/OufbXP32jZzL/+vCBxrWOLXaMSdbjuc2D+eZqm8OnRbVBZUQYH4Yc/H+Plne+we2w79PXjtRzl\n3uaHefH+p3hm3XW8U3NyMrh1uW/fXDj74INgzY7Nm+e2jRuhru6m/r58axl1HM46zmwYGwz30TbH\nWnpSKXrSabrTaXoiQawnnaYnDGPXerl2ORj5BR8/7+MXfLy8F9QXbLYYnlf05/bFyDEnqM8rO3be\n3nfCcnTv2nlttmQxCYNJBlssFQvKKUMsGZttN8kF9USkLfL52XIicizSNm+LR9rjc+3EufK8+IJ6\n+bzosfiCtqXUY/rBKLKcKazJsjY1BS//8jTfe7efnSPbudzST6zpLHfnHufFDVv5/Y1LeFLTWjh7\nNghk5XA2MACjo0EQ27w5uIW5aVOwjln6Bp76vAbX9xl2HAbDADYYCWLl8ojj0JRI0JNOszaWos9N\n0ltK0lVK0OEmaHVjrHLiZAvgz/h4016wn/Hwp8P9TGSfj9Tz/uzm5YO6iRli2RixTGxuXy6nI/VM\npJ6OYdImqIdtJhWpp+fqs+2psJwKP5uM1MthrFxWWBGRZUphTZaVQsHy/TeP8f2d23lvpJ/xuu0k\nsnnWp5/id+/dytce38qmro1Xf1LTdYOJ/wMDwRIZAwPBFo8HYeyhh+b2d94ZtH8Cju8zVCxyOl9g\naDzPufE8Y+eLTJ4vcnnSYXrSwbvk0VGM016MszofoyUfoyEPubwhPW1JTlvMtI9/2cOb8rC+JdGQ\nIF4fJ14fJ1YXI14XlON1wRarixHPheVcUJ7dZ2Pzy9kY8excOZaNEUtoEVgRkVtFYU1qWr7g8d23\n9vPD3dvZPdrPWK6fJDnuSm3lubuf4g+f3sqDPVd5UnNiIghk0e3oUVi7Nghj0a2z82P7Yq3Fu+zh\nnncpTZQojDuMj+YZHytwcazI1LhD8YKLO1mCSY/ERZ/sZUvDFGTyUMoavIYYNMVJNMVJNyfJNSWp\nb0mRbEoQb4yTaEwQb4gTb4gHgSwsz9br45iURplERJYThTWpKZNTBf7rzd1sG+hn3/ntnK/bQabY\ny13pp/id9Vv5xrNP8cDaNfM/5Lpw/HjwBGZ0u3gxuI354INz24YNs/PLrG9xJ1zccy7OqIM76uKM\nObhjLu6oy9RokemxIs64i3++RPyCh5cwTDXBxUaYqLc4zTFojpNYlSC7KkX96hTNrWlaWzN0tGXp\nbM+QXpUi0ZgI5h6JiIgsoLAmVe3D4Qv8x5s7eP3IOxye6udy3T5yM/dxT/YpXrh3K9/47JPc2R0+\nqWktDA0FT2MePBgEsoMHg6DW2xuEsY0bKd21EaftXpxYK85oieJwEWfEwRl2gv25YHPHS9iGGO7q\nODOrDJMtcL4Rhhs9Ttd7zLQYsm0pGtrTtHRkaG/PsKY5x5p0mjXhBP6U3hkpIiKfkMKaVI2S5/Pz\n3cf44c532XFmB6e8dylmT9My/Sgbm5/ixQ1P8vVnH6NzVX1wC/Pw4SCMHToEhw5hDx7Gja+iePuj\nFDo3UmxcRzHRjeM0Uhz1cIYcikNFrGtJdKcotcfJt8W5uNowvsoy1OTzUZPHiQaXMw0l6jvS9DRk\nWRuGr75MhrWZzGxda4eJiMitoLAmFXPk9Bjff2cXbx1/j8MXd3I+s5O420SPfZwtnZ/hiw8/zpce\n6CPzm+Nw+DDewFEKe85SODZJcaaewur7KOZuo2DbKE7XUZyIkWhKkOpN43clyHfGudBuGG2FMy0e\nJ5s9jjQ6HEsG64itzWToS6dnA1hfJjMbyjpSKb1LUkREqoLCmtwSg2OXePnXe/nVkT0cGN/NkHmP\nUmqclpkt3NPwKM933MPXGpvo+WCSwr5hCscvURgsUcg3U0ivpeC14nsp0p2QvD1HcV09l7tijHcY\nBlt9Tq72ONrsctIWGXVdulIp+jKZYDSsPCoWCWYaFRMRkVqhsCafuhOD53ll5wDbjw9wcHwPQ+zB\nzQzSdHkjz008xPNTa9g000TLGBTPeuQvZCl4bXjUkWmcxnTFKN6RY/LuJs7dUcepdsuJ1hJHckVO\nOQ6O78+Oii0MZH2ZDN2pFAnNFRMRkWVCYU1umON6vLX/JG8cOMjOU/s5fmmAKf991l3I8JkPH2LT\nxG2sm2qlaaoBphooOM0k41PQnKfYG2fyjiwj6xs4dVc9R3sTHGlwOO0G75vsi8wRWxjKVieTWppC\nRERWDIU1+Vi+b9l1bJA3DrzPex8c5vSZ42RHR7hn3OGekT76JvpYPdNOqrAa62cxqQmKrQ6Ta+Oc\nuyPLqTtyHF2XZX9PjJGET/ditygj88dyVfK+SRERkWqgsCazpvIObx84ybuHj/PhoQ8xJwbpGr7M\nmskEXZM9NM30kHTawYBTd5HJTp/RviSn12U4ekeG/X1J8m2G3rocayOT9aOBrDOVIq5RMRERkSVT\nWFthCk6JHYc+Ys/Ok4zt/ZDUiRFaRgqsupyjsdhG3UwnqWIDxVyeyVUOw72G07elObEuzZmeGPFu\nw+qOHGsaG2cD2RotZyEiInLTLMuwZox5AfgnIAZ821r7DwuOL9uwVvJ8Dh4b5kD/b7iw8yTeyTGS\nE5bcdJa6fBP1M83kphsoZEpMtJYY6jKcWptipBu8do9kT4yGvizdbatY09TEmkyG3jCItSQSmism\nIiJyiy27sGaMiQHHgc8CQ8Au4CVr7dHIOTUZ1opOiYN7Bjn1zjFGjgxRHJrGXPBJT6XJ5evITddR\nN11HupBgssVjpNMw3BnjcrNDoakIrZbG2zOsvruRrt52elpa6M5k6EmlaKqhIPb222/z9NNPV7ob\ncoP0/dUufXe1Td9f7VpqWKule1uPAiestacAjDHfA74AHL3mpyrk3OAkR3f+hjNHhrhwZpLCeAF7\n0ScxFSeVT5MuZMjOZMjNpGm4lKSYgcurU5Sa1pBvdMi3FCn2FWGVJdft03N/irX3d3FXVye/nc2y\nOplcdou76n84tU3fX+3Sd1fb9P0tf7UU1nqAM5H6IEGA+1RdnJji1MlzjJ4eZ2LkIlMT08xcyFO8\n5OBeLuHPWCgaYsUYiWKCpJMkVUyRLqbIFJJk8gly03FiPlxqBK8+RyKXJJF1KOSK5BsdSt0lYs0z\n5Dp82tal6d7Sw4b1a2jNZEhqHTERERGJqKWwtiT/et82wGDCu6HGGmK+IebHML4h7htiXoyYD/FS\njGTJkHQNSceQdIPPFDJQTFv8dIxkKkcumSGZLOGkXdxkCTfl4iVLFJpK5LMusfo86ZYUDW05sn2r\n6Livm/UP9FGXSVfsv4OIiIgsD7U0Z+0x4O+stS+E9b8EbPQhA2NMbfxmRERERGDZPWAQB44RPGAw\nDLwHfM1ae6SiHRMRERG5iWrmNqi11jPG/AnwC+aW7lBQExERkWWtZkbWRERERFaiZfPooTHmBWPM\nUWPMcWPMX1S6P7J0xphvG2POGWMOVLovcn2MMb3GmDeNMYeNMQeNMd+sdJ9k6YwxaWPMTmPMvvD7\n+9tK90mujzEmZozZa4z5SaX7ItfHGPORMWZ/+PfvvWueuxxG1payYK5UL2PMk8AU8B1r7cZK90eW\nzhjTCXRaaweMMfXAHuAL+rtXO4wxOWvtTDgv+NfAN6211/zBIdXDGPNnwMNAo7X285XujyydMeYD\n4GFr7YWPO3e5jKzNLphrrXWB8oK5UgOste8AH/uHVaqPtXbEWjsQlqeAIwRrIkqNsNbOhMU0wTzm\n2v8X/AphjOkFXgS+Vem+yA0xLDGHLZewttiCufqBIXILGWNuAx4Cdla2J3I9wtto+4AR4JfW2l2V\n7pMs2T8Cf44Cdq2ywC+NMbuMMX90rROXS1gTkQoKb4G+DPxpOMImNcJa61trNwG9wG8ZY+6rdJ/k\n4xljPgecC0e2TbhJbXnCWruZYHT0j8MpQYtaLmHtLLA2Uu8N20TkJjPGJAiC2n9aa7dVuj9yY6y1\nl4C3gBcq3RdZkieAz4fznr4LPGOM+U6F+yTXwVo7HO7HgB9zjVdoLpewtgu40xjTZ4xJAS8BejKm\ntuhfhrXr34D3rbX/XOmOyPUxxrQaY5rCchZ4DtDDITXAWvvX1tq11tp1BD/z3rTWfr3S/ZKlMcbk\nwjsSGGPqgOeBQ1c7f1mENWutB5QXzD0MfE8L5tYOY8x/AzuA9caY08aYb1S6T7I0xpgngD8Ang0f\nP99rjNHITO3oAt4yxgwQzDV83Vr7swr3SWQl6ADeCeeL/h/wqrX2F1c7eVks3SEiIiKyXC2LkTUR\nERGR5UphTURERKSKKayJiIiIVDGFNREREZEqprAmIiIiUsUU1kRERESqmMKaiIiISBVTWBMRERGp\nYgprIiIiIlVMYU1ERESkiimsiYiIiFSxRKU7ICJSrYwxfwXcD/w9sAVoAZqttX9T0Y6JyIqikTUR\nkUUYY54HtgFngH8B/gdoAn6vkv0SkZVHYU1EZHGrrLXvA48B/26tnQG+BXy5st0SkZXGWGsr3QcR\nkapkjMkCF4A7rbWDle6PiKxMGlkTEbm6J4HTCmoiUkkKayIiV/cs8GalOyEiK5vCmojI1d0N/KjS\nnRCRlU1z1kRERESqmEbWRERERKqYwpqIiIhIFVNYExEREaliCmsiIiIiVUxhTURERKSKKayJiIiI\nVDGFNREREZEqprAmIiIiUsUU1kRERESq2P8DFjTd8tx7AaIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1056d6438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(r, linear, label='linear')\n",
    "plt.plot(r, huber, label='huber')\n",
    "plt.plot(r, soft_l1, label='soft_l1')\n",
    "plt.plot(r, cauchy, label='cauchy')\n",
    "plt.plot(r, arctan, label='arctan')\n",
    "plt.xlabel(\"$r$\")\n",
    "plt.ylabel(r\"$\\rho(r^2)$\")\n",
    "plt.legend(loc='upper left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will show how robust loss functions work on a model example. We define the model function as \n",
    "\\begin{equation}\n",
    "f(t; A, \\sigma, \\omega) = A e^{-\\sigma t} \\sin (\\omega t)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which can model a observed displacement of a linear damped oscillator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define data generator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_data(t, A, sigma, omega, noise=0, n_outliers=0, random_state=0):\n",
    "    y = A * np.exp(-sigma * t) * np.sin(omega * t)\n",
    "    rnd = np.random.RandomState(random_state)\n",
    "    error = noise * rnd.randn(t.size)\n",
    "    outliers = rnd.randint(0, t.size, n_outliers)\n",
    "    error[outliers] *= 35\n",
    "    return y + error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define model parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = 2\n",
    "sigma = 0.1\n",
    "omega = 0.1 * 2 * np.pi\n",
    "x_true = np.array([A, sigma, omega])\n",
    "\n",
    "noise = 0.1\n",
    "\n",
    "t_min = 0\n",
    "t_max = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data for fitting the parameters will contain 3 outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_train = np.linspace(t_min, t_max, 30)\n",
    "y_train = generate_data(t_train, A, sigma, omega, noise=noise, n_outliers=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the function computing residuals for least-squares minimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(x, t, y):\n",
    "    return x[0] * np.exp(-x[1] * t) * np.sin(x[2] * t) - y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use all ones as the initial estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x0 = np.ones(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run standard least squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res_lsq = least_squares(fun, x0, args=(t_train, y_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run robust least squares with `loss='soft_l1'`, set `f_scale` to 0.1 which means that inlier residuals are approximately lower than 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res_robust = least_squares(fun, x0, loss='soft_l1', f_scale=0.1, args=(t_train, y_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define data to plot full curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_test = np.linspace(t_min, t_max, 300)\n",
    "y_test = generate_data(t_test, A, sigma, omega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute predictions with found parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_lsq = generate_data(t_test, *res_lsq.x)\n",
    "y_robust = generate_data(t_test, *res_robust.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x107031898>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAGECAYAAABtQ7cTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VFX+x/H3nfQEQkICIbSEFnqJkSYCkS4goohiWQs2\n1FXR1Z+uipRdXcsq2BZxpYkFsSwqNkQEbLSQ0HsJPRBSgCSkzf39cSEQasokc5N8Xs8zT5I7d858\nJ0TnM+fcc45hmiYiIiIiYi8OdxcgIiIiIudSSBMRERGxIYU0ERERERtSSBMRERGxIYU0ERERERtS\nSBMRERGxIVuGNMMwHjMMY51hGGsMw/jIMAxvd9ckIiIiUp5sF9IMw6gLPAxcZppmO8ATGOHeqkRE\nRETKl6e7C7gADyDAMAwn4A/sd3M9IiIiIuXKdj1ppmnuB14DdgP7gDTTNBe4tyoRERGR8mW7kGYY\nRhBwLRAB1AWqGYZxi3urEhERESlfdhzu7APsME0zBcAwjC+BK4CPzzzJMAxtOioiIiIVhmmaRnHO\nt11PGtYwZxfDMHwNwzCA3sDG851omqZu5XgbO3as22uoajf9zvU7rwo3/c71O68Kt5KwXUgzTXM5\n8DkQD6wGDOA9txYlIiIiUs7sONyJaZrjgfHurkNERETEXWzXkyb2FRsb6+4Sqhz9zsuffuflT7/z\n8qffecVglHSc1N0MwzArau0iIiJStRiGgVkJJg6IiIiIVHkKaSIiIiI2ZMuJAyIiInYSGRlJYmKi\nu8sQG4uIiGDXrl0ubVPXpImIiFzCyeuJ3F2G2Nil/kZ0TZqIiIhIJaGQJiIiImJDCmkiIiIiNqSQ\nJiIiImJDCmkiIiIiNqSQJiIiIgViY2Pp1atXsR6zePFixo/XltuuppAmIiIiBQyjWKtEALBo0SIm\nTJiA0+ksg4qqLi1mKyIi4gJOp5P4+HgAoqOjcThc0w9SVu260qn1wbSWnGvZ719aRESkgomPX09M\nzGh69EikR49EYmJGEx+/3rbtnjJ79mxatmyJr68vbdu2Ze7cuYXuz87O5vHHH6dt27ZUr16d8PBw\nhgwZwubNmwvOGT9+PBMmTADAy8sLh8OBh4dHwf3jxo0jJiaGGjVqUKtWLXr37s2yZctc9hoqM+04\nICIicgkXW03e6XQSEzOahIRJnO77cNKhw2ji4iaVuOerrNo9ZcGCBfTv359rrrmG+++/n8OHDzNm\nzBhyc3Np0aIFCxcu5OjRo/ztb3+jd+/e1K1bl9TUVP7zn/+wYsUKNm3aRO3atdm/fz9jx45l2rRp\n/P777wV1derUCYB7772XHj160LBhQzIyMvjwww/58ssviYuLo3Xr1qV6DXZSFjsOKKSJiIhcwsXe\ngOPi4ujRI5HMzOsLHff3/4IlSyKJiYkp0XOWVbundOvWjfT0dNatW1dwbNmyZXTt2pXY2FgWLlx4\nzmOcTifZ2dmEhYXxj3/8g0cffRQ43ZuWm5t70fDodDoxTZPWrVtz9dVXM3HixFK9BjvRtlAiIiJS\nak6nk5UrV3LDDTcUOt65c2ciIyMLHZszZw5dunQhODgYT09PAgICyMjIKDTkeTELFiygV69ehIaG\n4unpiZeXF1u3bi3y46syhbQqwOl0EhcXR1xcnGbeiIi4WHR0NFFRi4Az///qJCpqMdHR0bZrFyA5\nOZnc3FzCwsLOue/MY9988w0jRoygdevWfPLJJyxfvpyVK1cSGhrKiRMnLvk88fHxDBo0iMDAQKZN\nm8ayZctYuXIl7dq1K9LjqzrN7qzk4uPXM3LkFLZsiQUgKmom06bdT3R05bkOQETEnRwOB9Om3c/I\nkaPZsqUnAM2aLWLatFGlum6srNoFCA0NxcvLi6SkpHPuS0pKKuhNmz17Ns2aNWPq1KkF9+fl5ZGS\nklKk5/niiy/w8vLiyy+/LFRzamoqwcHBpXoNVYF60ioxp9PJyJFTSEiYRGbm9WRmXk9CwiRGjpyi\nHjUREReKjm5NXNwkliyJZMmSSFatesMlH4bLql2Hw0HHjh35/PPPCx1ftmwZu3btKvg5KysLT8/C\n/TkffPAB+fn5hY75+PgUnH+mzMzMQjM9ARYuXMju3btL+xKqBIW0Siw+Pv5kD9qZ/8wOtmzpWbDm\njoiIuIbD4SAmJoaYmBiXrmVWVu2OHz+eTZs2ce211/Ldd98xY8YMbrrpJsLDwwvOGTBgAJs2beLx\nxx9n4cKFvPzyy4wdO/acXrBWrVoB8O9//5vly5cTFxdX8Pjjx49zxx13sHDhQiZPnsxf/vIX6tev\n77LXUZkppImIiFRBvXv35qOPPmLLli0MGzaM1157jTfeeIPmzZsX7Dpw77338uyzzzJnzhyGDBnC\nDz/8wLx586hRo0ahnQkGDx7Mgw8+yOTJk7niiisKlt/o168fb775Jn/88QfXXHMNM2bMYNasWTRt\n2rREOxtUNVqCoxIr6zV2RESqikstryCiddLOoJBWNKcnDpy+6HT69FGaOCAiUgwKaXIpCmlnUEgr\nuoqw75uIiJ0ppMmlKKSdQSFNRETKi0KaXIp2HBARERGpIhTSRERERGxIIU1ERETEhhTSRERERGxI\nIU1ERETEhhTSRERERGxIIU1ERETEhhTSREREqpivvvqKiRMnursMuQSFNBERkSpm7ty5CmkVgEKa\niIiInFdOTo67S6jSFNJERESqkLvuuouZM2eyb98+HA4HDoeDxo0bs3jxYhwOB//73/+47777qF27\nNnXq1AHgzjvvpFGjRue0FRsbS69evQodS05OZtSoUdSvXx9fX19atmzJf//733J5bZWNp7sLEBER\nkfLz/PPPc/jwYVauXMk333yDaZr4+PiQlpYGwCOPPMLVV1/Nhx9+yIkTJwBr30nDOHfbybOPHTt2\njG7dupGdnc2ECROIjIzkxx9/5IEHHiAnJ4eHHnqo7F9gJaKQJiIiUoU0atSIWrVq4e3tTceOHQuO\nL168GIDOnTvz3nvvlajtSZMmsWfPHtatW0fjxo0B6NWrF6mpqYwfP54HHngAh0ODeEWlkCYiIuJC\nxvhze5xcxRxrllnbpwwdOrTEj/3xxx/p3LkzERER5OfnFxzv168fU6dOZcOGDbRp08YVZVYJCmki\nIiIuVB5BqiyFh4eX+LGHDh1i+/bteHl5nXOfYRgcOXKkNKVVOQppIiIiUuB81575+vqed6bnkSNH\nCA0NLfg5JCSEsLAw3nzzTUzz3LDavHlz1xZbySmkiYiIVDE+Pj5kZWWdc/x8AQ0gIiKCpKQkjhw5\nQkhICADbt29n8+bNhULagAEDePvtt2nQoEGh41IyunpPRESkimnVqhUpKSm8++67rFy5knXr1gGc\nt/cLYPjw4QDceuutzJ8/n48++oihQ4dSq1atQuc99thj1K5dmyuvvJIpU6awaNEivv32W1577bVS\nXetWVaknTUREpIq55557WLZsGc8++yxpaWlEREQwffr0C/akNWnShC+++ILnnnuO6667jqioKCZO\nnMiLL75Y6DGBgYH88ccfTJgwgVdeeYV9+/YRFBRE8+bNGTZsWHm9vErDuFBqtjvDMMyKWruIiFQs\nhmFcsJdJBC79N3Ly/mJN/dVwp4iIiIgNKaSJiIiI2JBCmoiIiIgNKaSJiIiI2JBCmoiIiIgNKaSJ\niIiI2JBCmoiIiIgNKaSJiIiI2JBCmoiIiIgNKaSJiIiI2JBCmoiISBUzbtw4HA5FALuz5b+QYRg1\nDMP4zDCMjYZhrDcMo7O7axIREaksDMO44GbqYh+e7i7gAt4AvjNNc7hhGJ6Av7sLEhERESlPtutJ\nMwwjEOhumuZ0ANM080zTPOrmskRERCqtN954g1atWuHv70/NmjXp2LEjX331VcH9TqeT5557jrp1\n6xIQEECvXr3YsGEDDoeDCRMmuLHyys2OPWmNgGTDMKYD7YGVwKOmaWa5tywREZHK56OPPuKJJ55g\n3LhxXHnllWRlZbFmzRpSUlIKzhk7diz/+te/eOKJJ+jbty8rV65kyJAhGjItY3YMaZ7AZcBDpmmu\nNAxjEvA0MNa9ZYmIiFQ+S5cupX379jz77LMFxwYMGFDwfVpaGpMmTWLUqFG8/PLLAPTp0weHw8HT\nTz9d7vVWJXYMaXuBPaZprjz58+fAU+c7cdy4cQXfx8bGEhsbW9a1iYiIXFxZ9i6Zpsub7NixI5Mn\nT+aRRx7h2muv5YorrsDPz6/g/rVr15KZmcnw4cMLPW7EiBEKaRexaNEiFi1aVKo2bBfSTNNMMgxj\nj2EYUaZpbgF6AxvOd+6ZIU1ERMQWyiBIlaXbb7+d7Oxspk6dyuTJk/H09GTgwIG8/vrrREREcODA\nAQDCwsIKPe7sn6WwszuPxo8fX+w2bDdx4KRHgI8Mw0jAui7tRTfXIyIiUmnde++9LF26lOTkZD74\n4AOWL1/OiBEjAAgPD8c0TZKSkgo95uyfxfVsGdJM01xtmmZH0zQ7mKZ5vWma6e6uSUREpLKrUaMG\nw4cP58Ybb2TdunUAtGvXjoCAAObMmVPo3E8++cQdJVYpthvuFBERkfJz//33U716dbp27Urt2rXZ\nvHkzs2bNon///oAV3B577DFefPFFqlWrRr9+/VixYgVTp07V7M4yppAmIiJSBZ0KWN26dWPGjBl8\n+OGHpKenU7duXW6//fZC132f+v7999/nnXfeoUuXLsybN49WrVq5ofKqwzAr2AWOpxiGYVbU2kVE\npGIxDAO955zL4XAwbtw4nn/+eXeX4naX+hs5eX+xuh5teU2aiIiISFWnkCYiIiIloo3ay5auSRMR\nEZESyc/Pd3cJlZp60kRERERsSCFNRERExIYU0kRERERsSCFNRERExIYU0kRERERsSCFNRERExIYU\n0kRERERsSCFNREREis3hcJTrdlCrV69m/PjxpKWlFen8yMhIRo4cWcZVlS2FNBEREbG9hIQExo8f\nT0pKSpHOrww7ISikiYiICDk5Oe4u4aJM06wUwas4FNJERESqmHHjxuFwOFi/fj0DBgygevXq3HTT\nTQX3T5w4kRYtWuDj40PdunV5+OGHOXbs2DntmKbJiy++SIMGDfD396dnz56sXr260DkXGnZ0OBxM\nmDCh4OetW7dy3XXXERYWhp+fHxEREdx00004nU5mzpxZ0EbTpk1xOBx4eHiwe/fuIr/mpKQk7rjj\nDurVq4evry9169ZlyJAhJCcnF5yzc+dOBg0aREBAAGFhYYwePZopU6bgcDiK9Vyuor07RUREqphT\nPVJDhw7l7rvv5umnn8bhsPptnnnmGV566SUefvhhBg8ezIYNG3juuedYs2YNixcvLtTOzJkziYiI\n4J133iE7O5sxY8bQp08ftm7dSlBQUKHnupSBAwcSEhLClClTCAkJYd++fXz33Xc4nU4GDRrEc889\nxwsvvMAXX3xBvXr1AAgPDy/ya77tttvYs2cPr732GvXr1ycpKYmff/6ZzMxMAHJzc+nTpw/Z2dlM\nnjyZWrVqMWXKFL788kv39eCZplkhb1bpIiIiZa+yveeMGzfOdDgc5ltvvVXoeEpKiunj42OOHDmy\n0PEPP/zQNAzD/OabbwqOGYZh1qpVy8zKyio4tmvXLtPLy8t8/vnnC45FRkaad9111zk1GIZhjh8/\n3jRN00xOTj6n/bPNmDHDdDgc5vbt24v0Gs9+3mrVqp3zes/03nvvmQ6Hw1y+fHnBMafTabZu3dp0\nOBxmYmLiRZ/vUn8jJ+8vVtZRT5qIiIgLGYsWlVnbZmysS9sbOnRooZ+XLl1Kbm4ut956a6HjI0aM\n4K677mLx4sUMHjy44PjAgQPx9fUt+DkiIoIuXbrw559/FquOkJAQGjduzNNPP83BgweJjY2ladOm\nJXhFF9axY0deffVVnE4nvXr1ok2bNoXuX7p0KQ0aNKBjx44FxwzD4MYbb2T8+PEuraWoFNJERERc\nyNVBqiydPVx4aubk2cc9PDwICQk5Z2ZlWFjYOW2GhYWxYcOGYteyYMECxo0bxzPPPENycjKNGjXi\nySefZNSoUcVu63zmzJnD+PHjefXVV3nssceoU6cOo0aNYsyYMQAcOHDggq/HXTRxQEREpIo6+1qr\nmjVrYpomBw8eLHQ8Pz+fI0eOULNmzULHk5KSzmkzKSmp4JoxAF9f33Nmjp5vGY3IyEhmzJjBoUOH\nSEhIoHfv3jz44IP8+OOPxX5d5xMaGspbb73Fnj172LRpE3fddRdjx45lypQpgBVMz/d6zv5dlCeF\nNBEREQGgS5cueHt7M3v27ELHZ8+eTX5+PrFn9RJ+9913ZGVlFfy8a9culi5dyhVXXFFwLCIignXr\n1hV63Lx58y5aR7t27XjttdcACh7r4+MDUOj5SqpZs2b885//JDg4uKD9rl27smfPHpYvX15wnmma\nzJkzp9TPV1Ia7hQREREAgoOD+dvf/sZLL72Ev78/AwcOZMOGDYwZM4bu3bszaNCgQuf7+fnRr18/\nnnjiCU6cOMHYsWMJCgpi9OjRBeeMGDGCu+++m8cff5zBgwezevVqZsyYUaidtWvX8uijj3LTTTfR\ntGlT8vPzmT59Ol5eXvTq1QuAVq1aYZomb7/9NnfccQdeXl60b98eT89LR5mjR4/Sp08fbr31Vlq0\naIGXlxdz584lLS2N/v37A3DHHXfw0ksvcf311/PCCy9Qu3Zt3n33XY4fP17K32rJKaSJiIhUQRda\nVuLMgDJ58mRCQkK48847efHFF895/B133IG/vz9//etfOXLkCJ06deKzzz4rWH4DrPCzd+9epk6d\nynvvvUePHj2YO3cuTZs2LaihTp06REREMHHiRPbu3Yuvry9t27bl22+/JTo6GrB618aPH897773H\n+++/j9PpZOfOnTRs2PCCr+9U+76+vsTExPD++++TmJiIw+GgefPmfPzxxwUTIby8vFiwYAF//etf\neeihhwgICOCWW25h0KBBPPDAA6X7ZZeQYc0KrXgMwzArau0iIlKxGIaB3nOqplML6V4sEMKl/0ZO\n3l+sBdd0TZqIiIiIDSmkiYiIiNiQhjtFREQuQcOdcika7hQRERGpIhTSRERERGxIIU1ERETEhhTS\nRERERGxIIU1ERETEhrTjgIiIyCVERERccIV+EbD+RlxNS3CIiIiIlDEtwSEiIiJSSSikiYiIiNiQ\nQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiI\niNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmk\niYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQbUOaYRgOwzBWGYbxtbtrERERESlv\ntg1pwKPABncXISIiIuIOtgxphmHUBwYC77u7FhERERF3sGVIAyYCTwKmuwsRERERcQfbhTTDMAYB\nSaZpJgDGyZuIiIhIleLp7gLOoxswxDCMgYAfUN0wjA9M07z97BPHjRtX8H1sbCyxsbHlVaOIiIjI\nBS1atIhFixaVqg3DNO07omgYRk/gb6ZpDjnPfaadaxcRERE5xTAMTNMs1uig7YY7RURERMTmPWkX\no540ERERqSjUkyYiIiJSSSikiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqI\niIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQ\nQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiI\niNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDSmkiYiIiNiQQpqIiIiIDXkW5STDMOYBB4GFwELTNA+W\naVUiIiIiVVxRe9JeBzKA/wP2GoaxwTCMtw3DuN4wDP+yK09ERESkajJM0yzeAwwjGOgBjAAGA7nA\ng6ZpznZ9eRetwyxu7SIiIiLuYBgGpmkaxXpMaYKOYRiPAr8DbwETTNP8vsSNFf+5FdJERESkQihJ\nSCvScKdhGOMMw1h18mvjM+5ymqa5EqtnrU9xnlhERERELqyo16R5AE8BkcBawzC2GIaxEuh68v7G\nwDbXlyciIiJSNRVpdidwADBN07zTMIy/At0Af+A7wzBqAOuAKWVUo4iIiEiVU+Rr0gzD6I4V1H47\nz31RwAHTNI+5uL6L1aNr0kRERKRCKPeJA+6kkCYiIiIVRZlNHBARERGR8lXUa9KkuNLSYMUKWL8e\nUlLA6YQ6dSAqCrp2herV3V2hiIiI2Jh60lzJNGHePBg8GBo2hBdegG3bwMMDvL2twPbii1C3LvTt\nC599Brm57q5aREREbEjXpLnK0qXw8MOQnQ1PPgnXXguBgec/9/hx+O47ePtt2L3bCm4jRoBDmVlE\nRKQy0sQBd8jNheeeg5kz4fXXix+2fv0VHn/c6mmbOROaNi27WkVERMQtNHGgvKWmwoABsG4drFkD\nt9xS/N6w7t1h2TK46Sbo0gWmTi2bWkVERKRCUU9aSR04AL16wdVXw6uvWtedldbGjXDdddCnD0yc\nCF5epW9TRERE3E7DneXl0CGIjYXbboNnnnFt2+npcOutkJcHX3wBAQGubV9ERETKnYY7y8OxY9Cv\nH9x4o+sDGkCNGjB3LoSHWzNAU1Nd/xwiIiJie+pJKw6n0xqOrF0b3nsPjGIF4uI/12OPWbNGf/rp\nwjNFRURExPbUk1bWnn/e6tl6552yDWhgTUCYNAkuuwwGDYLMzLJ9PhEREbEV9aQV1c8/w+23Q3y8\n1ZNWXpxOuPNO61q1L790zQQFERERKVfqSSsryclwxx0wfXr5BjSwetTef9/qSXv4YWtXAxEREan0\nFNKK4qGHrIkC/fq55/m9veHzz2HJEpgyxT01iIiISLnScOelzJsHo0fD2rXg51f2z3cxW7dCt27w\n9dfWwrciIiJSIWi409WOH7d60aZMcX9AA2jWzBr6HD7cWqtNREREKi31pF3MU0/BwYPWnpp28txz\n8Pvv1tIcnp7urkZEREQuQTsOuNKOHdCxo7UvZ3h42T1PSeTnW8tytGsHr7zi7mpERETkEhTSXOmG\nG6BDB6vXyo6OHIH27a1evt693V2NiIiIXIRCmqv89hvccgts3myPa9Eu5Mcf4d57YfVqCA52dzUi\nIiJyAQpprmCacNVV1rpod93l+vZd7eGHISUFPvrI3ZWIiIjIBWh2pyv88gvs3w9/+Yu7Kymal1+G\nVatg9mx3VyIiIiIuZLuQZhhGfcMwFhqGsd4wjLWGYTxSbk9umjBmDIwdW3FmTfr7w6xZ8MgjsHev\nu6sRERERF7FdSAPygMdN02wNdAUeMgyjRbk8888/W0OHI0aUy9O5zOWXW8Oe992nbaNEREQqCduF\nNNM0D5qmmXDy++PARqBeuTz5yy9ba6NVxE3Mn3oK9uyBTz91dyUiIiLiAraeOGAYRiSwCGhzMrCd\neZ9rJw7ExcHQobB9u7VXZkW0bJn1Gtatg5AQd1cjIiIiJ1WqiQOGYVQDPgcePTuglYlXXoHHH6+4\nAQ2gc2dry6gnnnB3JSIiIlJKtrw63jAMT6yANss0za8udN64ceMKvo+NjSU2NrZkT7hzp3U92tSp\nJXu8nbzwArRubb0eLXIrIiLiFosWLWLRokWlasOWw52GYXwAJJum+fhFznHdcOeTT1pfX33VNe25\n27x5MHo0rF1r78V4RUREqohKsZitYRjdgCXAWsA8eXvGNM0fzjrPNSEtIwMiImDFCmjUqPTt2cWN\nN0KLFjBhgrsrERERqfIqRUgrKpeFtClT4PvvYe7c0rdlJ3v3WnuPLlsGTZq4uxoREZEqrVJNHCgX\npglvvWWtMVbZ1K9vDeM++qi7KxEREZESqNoh7Y8/IDcXevVydyVl47HHYOtW+OYbd1ciIiIixVS1\nQ9p//wv33ANGsXofKw5vb6un8NFHISvL3dWIiIhIMVTda9LS0iAyErZsgdq1S12P0+kkPj4egOjo\naBwOG+XfYcOgXTtrT1IREREpd7omrTg++gj69XNJQIuPX09MzGh69EikR49EYmJGEx+/3gVFusjE\nifDmm5CY6O5KREREpIiqbk9adLS1V2e/fqWqw+l0EhMzmoSESeBlQvNj0OQ4Ya3ncdOdvfF2OKjj\n7U0zPz86BwYS5q4dDcaOtXoNP/nEPc8vIiJShWkJjqJaswYGDYJdu0q9mfqylSvp8fcD5MRGwmWp\nkBgAW6vhlbadv94TTFi9ehzIzmZzVhZ/pqfT0NeXm2vX5i9hYdT39b1guy4fPs3IgObN4bPPoGvX\n0rUlIiIixaKQVlRPPAFeXvCvf5X4+Z2mycyDB3l+yxb2r83BOa8N/FoLMq2dtvz9v2DJkkhiYmIK\nHpPndLLs2DFmHTzIZ4cPMzgkhGcjIojy9y/Udnz8ekaOnMKWLbEAREUtYtq0+4mObl3iegGYORMm\nT7ZmtdrpmjkREZFKTiGtKPLyoEEDWLgQWrYs0XOvOX6c+7dswQG8EBnJY73HWcOdBZf4OenQYTRx\ncZMu2AOWlpvLW/v28ea+fdweFsbYyEgCPT0LD58Wo70icTqhUydrI/lbbil5OyIiIlIsCmlF8cMP\n8PzzsHx5sR9qmiZv79vHPxITeaFRI+4OD8dhGAU9X5u39MAMPEx4+2+5+f4YgsOq42F4EOofSmRQ\nJO3rtCfQJ7BQm4dycnh6xw4WpKYyo0ULauzYQY8eiWRmXl/ovPP1zJXIr7/CbbfBpk3a11NERKSc\nlCSkeZZVMbY1axbcfnuxH3YiP5+7Nm9mS2Ymf152GU1OBhzTNEkPPkzHCSfYu/4BTEyi6l/OYc8D\nnDh2nDxnHiv2r2Bn2k7WJq2leWhzhjYfyog2I2gW0oza3t5Ma9GCH1NS+MvGjfQBzLIciezeHTp2\nhNdfh2efLcMnEhERkdKoWj1pWVkQHg6bN0NYWJEflpKby9B16wj39mZmixb4enjgNJ18svYTXvr9\nJQBub3dmqHLVAAAgAElEQVQ717W8jibBTTAusDhudl42y/Yt4/MNnzNn/Ryiw6N58oon6dXI2vEg\nOSeHERs2sGzZZo7/31A46nPykS4a7jxlxw5r2HPtWuv3ISIiImVKw52X8r//WSvwL1xY5Iccyc2l\nd0ICVwUH81qTJjgMg2V7l/HX7/+Kp8OTcT3H0a9JvwsGsws5kXeCj9d+zL9++xeNghrxat9XaV+n\nPXlOJ3ctXcFnew5hPJ2H45CTZs0WMX36qNJPHDjT//0fHDkCU6e6rk0RERE5L4W0S7nlFujRA0aN\nKtLpqbm59F69mj7BwbzcuDF5zjzGLRrHtIRpvNr3VW5te2uxw9nZcvNzmRI3hQmLJzAyeiRje47F\nz8uPN/fs4Z87dvCWnx/DL7/c9TsYpKdbS3J8/721ZpyIiIiUGYW0izk11FnEbaBO5OfTd80aLqtW\njUlNm3Lw+EGun3M9Nf1qMm3INMKqFX24tCgOHj/Iw98/zPpD65kzfA5tardhdlISj2/fzoL27WkV\nEODS5wPg3Xfh00+tnsXKun+piIiIDWhbqIv5/nuIiSlSQHOaJnds2kRdb28mNm1KwsEEOr/fmYFN\nBzLv5nkuD2gAdarVYc4Nc3iq21NcNfMqZiTMYERYGK82aULf1avZlJHh8ufknnsgKQm+/db1bYuI\niEipVJ3ZnZ99BjfeWKRTx+/axb7sbBa0b8/vu39j2Jxh/GfQf7ih1Q1lWqJhGNzR4Q461evEkNlD\nWJO0hlf6vkKeadJ3zRoWtm9Ps7MWvi0VT0949VVrcd8BA6yfRURExBaqxnBnMYY6v0lO5sGtW1kZ\nE8O6vb9y8xc38/Gwj+nTuI8Lqi661KxUhn82nCDfID66/iM+OpzCuF27+D06mgYX2U6q2EwT+vSB\n4cOLfK2eiIiIFI+uSbuQL7+E//wHFiy46Gk7srLosmoVX7VpA0c3MGT2EL648Qt6RPRwQcXFl52X\nzS1f3kJGTgZf3vQl7xxI5sOkJH6Ljqa6K3u94uNh4EBraZLAwEufLyIiIsWia9IuZM6cSw515jqd\n3LJhA39v2JDq2XsY+ulQZl03y20BDcDH04dPb/iUsGph9P+wP/fWCqRzYCAjNmwgz+l03RNFR0O/\nfvDKK65rU0REREql8vekZWZC3bqwdSvUqnXB08bs3MmKo0eZ3jiMrlO78EKvF7i13a0urLjknKaT\nR79/lD/2/sH3t/7Ebdv20Nzfn7eaNXPdk+zZAx06wOrVUL++69oVERER9aSd1w8/WNsgnSegOZ1O\n4uLimLF8Of/dv593mzXiuk+HcmeHO20T0AAchoM3r36T7g27c/2nQ5gR1YiFqam8uXev656kQQPr\nmrTnnnNdmyIiIlJilb8n7c47rZD20EOFDhdsir4zlhNv1aDhwpVc1nMlngEGn97waakXqS0LTtPJ\nHXPvICUrhYnXzqZ7who+a92aHkFBrnmCo0etBW6/+04L3LqQ0+kkPj4egOjoaNcvTCwiIraniQNn\ny8+3ZnWuWAEREQWHnU4nMTGjSUiYBPfuhHpZMHc5Pr2fIOkfO6nhV6OMqy+53Pxcrvv0OoL9grm5\nx0Tu3byFuJgY6vj4XPrBRfHuu9ZyJQsWaIFbFzj1YWDLllgAoqIWMW3a/a7d4ktERGxPw51nW74c\n6tQpFNAA4uPjrTfNyEwYeBA+NKHvUxhznmfbhm3uqbWIvDy8mDN8DrvSdvHTihe4Jzycm1w5keCe\ne2D/fqs3TUrF6XQycuQUEhImkZl5PZmZ15OQMImRI6fgdOXEDylXpy6TiIuL07+jiJSpyh3SvvkG\nBg8+710mwMPbYFY96HMbLHgJx5EG5VpeSfl7+fPNzd/ww/YfqJM8H1+Hg2d37nRN46cWuH3yScjL\nc02bVVTBh4FC/5k52LKlZ8Hwp1Qs8fHriYkZTY8eifTokUhMzGji49e7uywRqaQqd0ibNw+uueac\nw9HR0dS+cSvUyIWMtyGlKcTfSVTUYqIryLVYQb5BfD3ia8Ytfp4HAlKYfegQcw8fdk3jgwZZPZBT\np7qmPZFKQD2jIlLeKm9IS0yEgwehU6dz7soyTbLv7Er932ZhtJ2G34JraN9+NNOm3V+hLupuFtKM\nT4Z9wqi5N/Pv+oHcv2UL+7KzS9+wYcC//w3jxsGxY6Vvr4qKjo4mKmoRcOYbuLNCfRiQ09QzKiLl\nreIkkuL65htrFX0Pj3PuejExkZ6hNQm8/Gf+1eMpfp3fllWr3qiQF3P3atSLsT3HMvbr4dwTFsrt\nGzfidMVkkMsug759tcBtKTgcDqZNu58OHUbj7/8F/v5f0L79oxXuw4CIiLhH5Z3dOWAA3HsvDBtW\n6PC2zEy6rFrFnXm/seXAH3w14itbLrdRXA9++yB7ju4jrflYrgkN5f8aNix9o7t3W0txaIHbUtES\nHJVDoVnhBZ9vnXToMJq4uEn6dxWRi9ISHKccOwb16sHevefsRTl07VqaeGTzwbdXE39/PPUDK0f4\nyMnPoeeMnvSMupFpjk5817Ytl7tiH85nnoEDB2D69NK3JVLBnV5SpScAzZotYvr0URWyF15EypdC\n2in/+x9Mngzz5xc6vDQ9neEb1tNg0zPc2no4D3V66PyPr6D2pO+h4387ct+A2czO8GNVTAzVSrsR\n+9GjEBVl7dzQoYNrChWpwNQzKiIloXXSTrnA0hvP7txJTxLBmcMDHR9wQ2Flq0GNBsy6bhbv/3gL\n7f28eGrHjtI3GhgIzz8PTzwBFTTQi7iSw+EgJiaGmJgYBTQRKVOV7/8wTid8++05IW1BSgqJJzL5\n8ffHmTJ4Cg6j8r10gL5N+jLq8lHsjX+ar5OT+Tk1tfSN3nsv7Ntn9aaJiIhIuah8SSU+HmrWhMaN\nCw6ZpskzO3cSdfQPbmgxlLZhbd1YYNl7rsdz1PBw0O3EMu7etImjpV2U1svLmuX5xBNa4FZERKSc\nVL6QNn8+9OtX6NDc5GSO5WaxfPUrTLhqgpsKKz8Ow8HMoTNZEv8qrb2yeWL79tI3Ongw1K6tCQQi\nIiLlpPKFtJ9+KhTS8k2T53buJHD/Zzx1xZPUCqjlxuLKT1i1MKZfO53Vv9/HD0eS+TElpXQNnlrg\nduxYLXArIiJSDipXSMvIgBUroGfPgkOzDx3CyM/g8N55PNL5ETcWV/76N+3PTS2uoUHSp9yzeTNp\nubmlazAmBnr3tvb2FLew8+bedq5NRKQiqlwhbckSK0hUqwaA0zR5MTGRjK2T+XffV/Hx9HFzgeXv\nxd4vknXoNxrlH+QxVwx7vvACvPOONZFAylVZbO7tqmCljcdFRFyvcoW0+fOtrYxO+jo5mYwTKTRw\nJnFdi+vcWJj7+Hj68MmwT9iw7BEWHDnEvOTk0jXYsCHcdx+MGeOaAqVIymJzb1cFK208LiJSNipX\nSDvjejTTNHkhcRcZ29/j5d4vVYqtn0qqeWhzXuk1Hq9tE7lvy2ZSSzvs+fTT8N131nZRUi5cvbm3\nK4OVNh4XESkblSek7dtnbV902WUA/Jyayp6MFDr55NC1QVc3F+d+d3W4i44B3oRmbOTx0g571qhh\n9aRpgdsKS8FKRMT+SrlnkI0sWGBd1O7hAcA/EneSvXMa/+hf+ZfcKArDMJgyeArt3uvM976N+aFW\nLQaEhJS8wfvugzffhB9/tDazlzIVHR1NVNRMEhKGcubm3lFRi4mOLtuh/HzTZGdWFtvOuB3OzeV4\nfj5peTkcy8+Ft/LgSDyc8IIsD0jyJsQ7kcMRPdmVlUWEr2+V7s0WESmJyhPSzrge7c/0dNamJ3OV\nv8ll4Ze5uTD7CPINYtaQ97h+/gvc4/E8Gzp1IbCke3ueucBtnz5Q2j1C5aIcDgfTpt3PyJGjC23u\nPW3aqBJtTXSx0Nem/bUsTkvjt/R0fklJZunRo3iZOQTkpUDWfrKObSMzYx/ZOWl4mrn4OsCZ7wRv\nE3z8ICgYzwbh5IRHce/qhRz1CMZweBFTPZAuNYLoERRE9xo18D/5gUpERM6vcmyw7nRCeDgsWwaR\nkQxcHc+vCa+z9OqnaV27tXsLtaG/L/g7H+eEM6BJf6Y0b17yhkwTrroKbrsN7rnHdQXKBblyc+/4\n+PWMHDmFLVt6YnpCnUGbaPVgT5bk5+CTl4ozNZ7M5GVcFuBN17CWNAluQpOaTWgU1IiwamEEeAXg\n4TgdtPLz81mxagXHc48TFBHE3mN7SUxLZGPyRpYd3saGE/lUD4nBs2ZH0j1r0TGwOteEhnFdaChN\n/f3L7HWKiNhBSTZYrxwhLSEBbrwRtmxhfUYGXVb8zqCUj5h9/Uz3FmlTOfk5dJp+FbubjePzdpfT\nKzi45I2tXAlDhsCWLQVLn0jFsSUjg7HxK/kqJwe/3CQy9v9AR58shjfpSfeG3Wkb1hZPh2t6SXPz\nc9mYvJHFuxbz/a4lLE5Lw7d2T7KDO9HQ14/bwxsyonZtUjftOBkeYwGIilrEtGn3Ex2tD1wiUnFV\n3ZD26quQmAhvv80dG9bwefxkEgY8TrOQZu4t0sY2JW+i85ePENh6DBs7d6VaaYYrb7sNmjaFceNc\nVp+UHdM0+Sk1lQnbNxB3PAPzwHdc4ZHM3S0HMrDZQIL9ShHaiyHPmceyvcv4fOOXfLxnPdk1ryAn\nuAuO7YfJ+KgP/F4Lch2Akw4dRhMXN0k9aiJSYVXdkNavHzz0EIeuvpqIP5ZwTfpnzLl2insLrADe\nXfkuz+47xs0th/B2VCmGPRMTrVm1a9ZAvXquK1BcyjRNfkpJYfTm1ew6fgS/A//jgYiW3B99Jw1q\nNHB7bSv3r+TfP7/JnCOpUGcYVGsI39aHzyPxz5rHkiWRxMTEuLVOEZGSqpohLTsbQkNh716eS07i\ntVWzWBV7Cy1rtXR3ibZnmiYDPr2RP8Lu5vvoLlwZFFTyxp59FnbuhI8/dl2B4jJ/pqdz97oV7MxI\nIeTQd7zcrj83tRnusqFMV4mLi6N77DayGnhC7JfQqR6E9cWxLJOPB4ZzU6dO7i5RRKREqmZIW7IE\nnniC7KVLqbXkZ7oc+Yz5w95zd3kVxqGMQ7SYfQ/Vmj/O5q7d8SvpjLuMDGjVCmbOhNhYl9YoJXcw\nO5u71i7j57RUwg9/y2vtr+b6ltfiMOw5bOh0OomJGU1CwiTAATW3QY/JGAOT8WlyIx2qBfKPpq3p\nHRysJT3EVjTZRS6lJCHNXh+jS+KXX+Cqq/jgwD5yjm7kpStGubuiCqV2QG0+7DGK4euW8tS22rzZ\nvFXJGgoIgNdfh4cfhlWrrCU6xG3ynE7Gbl3La/v24314Ae+0uIy7e79r23B2yjlLjZyAZjtzeLP9\nKJbmLeWNtXEMOzyYBp4B/PuEk/5JSRipqXDixOlbfr7193fmzdvbWoQ5JOT0rWZNqFUL/Pzc/bLF\nTVwVrE7PlI4FICpqpia7iEtU/J60nj0x//53GvobhB36ipU3/MfdpVVII799jE/8+rL48ivpFBhY\nskZME/r3h4EDYfRo1xYoRZZwNI3Bq37lUPoOHqiRy7+6PYi/l/+lH2gXubk4N2xg17ff4rdjB3VS\nUjB27IDERMycHNLqBDO5SzumDLmZWk54MXEffTMzMfz8wNcXHA7IzS18y8mBtDQ4cgRSUqyvR47A\n4cNWeIuIKHxr1gxatrT2qlWPSKV0brAq2Szic3p/raO2muyiXj57qHrDnZmZUKsWP23eyKC1v7Kg\nRX16RPZwd2kVUmZuJs1m348ReTvbu/XGp6T/EW/aBN27w9q1UKeOS2uUi8tzOvnruj94/1AqrY7/\nztc9RxEZHOnusi4tKQkWL4alS621DhMSoEEDaNPGGkJv3doKTQ0bWj1ghoHTdDJ30zc8kTCPfTX7\n0DggmLdaxdCn5rm7aFz0DcrptJ4/MfH0bdcu2LoVNm6E1FRo3twKbK1aQXS0NUlGf9sVmiuDVVxc\nHD16JJKZeX2h4/7+X9hisourwqiUXpULafnz5+MYO5bLXh3DsQPz2XrDRHeXVaGt2r+Krst+5v4W\nfXmzZYeSN/TUU9Y+qh984Lri5KLWH02l98pFpB7bzWuR4TzUbrh9r9nKzrYuU5g/H37+2QpGPXpA\n167QuTNcfjkUsTfXNE3mbf2ex+K/YXdwL1pXq8G0tl2Jrl4dcMEb1NGj1gePjRth/XqIj7eG8318\nrLB25q1BA7Dr71wKcWWwsnNIqwi9fFVJlbsmbepfnqHNiO6sycrj83b93V1OhXdZ3ct4stZiXtm/\nl9vrNeLywBola2jMGKvX4eT1glK2Xtq8jOf2HKJ99nrW932AEP9S7MlaVo4ehW+/hblzrf1e27SB\nq6+GKVOsUFbCdfoMw+CaqIEMbnY187b9yIOrvqFTZjY9Av15r3UXRo6cUugNKiFhKCNHFuMNKjAQ\nOnWybqeYJuzebYW1Vavgv/+FuDjr+BVXWGHziisgJkbXu1UB7txX91Li4+NPfkA582/dwZYtPYmP\nj3d7L19VcWZvfnFV6J60X+nGfeP6kdwml6TrJ9i356ACyXfm0+bLv5Ee2pfEHlfjVdJPWl9/be3r\nuWaNdZ2QuNyx3Bx6/zGXVVl5TAjz5JnoG91dUmG5uVYg+/BD+P57uPJKuO46uOYaCAsrk6d0mk4+\nXPclj2/4k7SQXvBLOvlv3wBHvQvOKZMejlPB7c8/rdsff8CGDVYYPRXauna1etvE7TKyM4ge/Bhb\nM2+DuqlQ5ygEZuNfaweRLRuSaUKW6SAbBybGyRtw8qsBeOPE23DiY5iQncPB3cnkJPtDOtQ093HX\n4K50b96aljXCaFwt1C29Vnbu5asqzuzNz8wcVrWGO5N8ahL2vxk8lruR14f8n7tLqjQS0xKJWvQF\nI5t0YXLbK0re0LBh1vVEEya4rjgBYFnyHvqs+gP/nIMs7DqE1iGN3F3Sadu2wbvvwqxZ1k4Ut91m\nbdsWUn49fHnOPP7v638xMeU41I2FucEwsyOc8Ci/N6jMTGvbtFOh7c8/rWHSbt1O39q1K3Evolzc\n0ezjLD64id+P7Gb10RS2n8jhoNOLDI8aOL2CMJzZkJGKefQERvoJ/PJS6NAsnAY1axLk5U1NLx+C\nvX3xcjhwYOAwjIKvuc48juZmk56Xw9HcbI7n5XE0L4fdaamkmwbZ3r5k4MkJfMjxrA6GJ955aVQ3\nswh1OKnn40XzgEC61KxHrzrNqO9fijUqL0LDnSXniskW5/7+q9g1aU9fcy8v3RLF0ibd6dyxs7tL\nqlTeWT2HRw/78eflXekYFFqyRvbtgw4drIvCW5VwaQ85x6sbF/P03jT6ex7k6x534+lhgzf5/Hz4\n5hv4z3+sC//vvBPuu88KaW7idDqJ7vgwa+q0gVtSoEYH+Cic9ltmsmrFxPJ/gzJNK8D+/vvp2969\n1lDqqdDWpUuRr8eT0w5lpPDV3gQWHEokPiOTPU5fTvjUxct5giDzOA08oWVAADFBYXQLjaRtUB38\nPDzKbdbj/sxUVh5JZHXaQTYdS2XniUz25uRz2PTjhHcIDmcu1fPTCffIo5mfH5cFhdI/rBkeuw7j\n4XC4aHmQngA0a7aI6dNHaeLARbhqssW5PZmVJKQZhjEAOBU9p5qm+fJ5zjHDZ35AdsJCDv97qj4R\nlIEu3/6THb4t2X/VdXiW9Pf7zjswe7YV1PRvVCr5znwG/vYRC05U4+V6ATzR2gbXYWZkwIwZ1hp5\ntWpZ6+QNG2abIe5T/7PdvLMTuYOWk3d9JAGBjXilWXNGRbTE4e5LJI4csXrYToW2VausYHvllaeD\nW8OG7q3RZkzTJCF5G7N2reLnlMNsyfPjhG89/PKP0dCRRXS1avStHcGQeq0I9bH/NYH5znxWHtnJ\n4qQdxKUfYnNmBjtPmBz1CgHPapByEN8jB+gVWYfro9oyuF4rwnwDivUcWoKj6Mp25m8lCGmGYTiA\nLUBvYD+wAhhhmuams84z/ed+yHd1mtCzcxc3VFr5pWWlEb7gY0bUbcz0mAElayQ/33qjGTnS6lmR\nEknNzqDDLx9zBB9+ubw7HUPdPLyZlgaTJlk9Z926wZNPWtdd2dCZb1ARLSJ4eOkMPjsRRIhfCG82\nb8dN9Zq4ucIz5ORYQe1UaPvtNw2RAuuO7ODdbUv5MSWZXWYg+T5h1MpP4TJ/L4aGN+amBm0J8vZx\nd5kucTokTIQ6+yF6A7Tcj3eLvXjUq0+WTx288rOoxTFa+HrSNbg2g+u2oFPN+u7/0FEJuPI6vko5\n3GkYRhdgrGmaV5/8+WnAPLs3zTAMs9e3b/DzwEfcUWaVMWf7r4zYmcavHdrTLbSEn+jXroVevaw3\nH104XWzrU/fRZelP1HTksfqqWwjyduPCtOnp8MYb8Oab1gSAv/8doqLcV08JHTh2kHv+/IAf8uvT\nwMeb/7bpTN/aNvzbrKJDpKkn0nl/6xI+P7CLNTneZPs1JMyZwpXV/bi5fisGhTfDp6Rb2NncpUJC\nq3at+Wn/OuYn7WDl0RS25UCKR03wrE6NvBQaeTm5PDCIvmGNGRAeRXXPIuz+cmof7FO33FxrDcH8\n/Et/NU1rlOTUzTAK/3yhm7f36d1AvL2tDx42CJmunmxx5nBzZuYNlSKkDQP6m6Z538mfbwM6mab5\nyFnnmQnJO2hvpwumK6n+C6ewIq8ah/rcXPJhz3/+0+oV+P57W/yHWFF8vWc1wzZspqt3Jr/0vB0P\ndw1THDsGb70FEydaO0qMGePW681cZUfabu5Y+gm/O5rS0tvJBx16EhNU291lXdz5hkibNYOOHa1l\nP2JioG1bqweuAtmRvpfXNv7CV8lH2O8dib+Ry2U++Yyo25TbI9tTrShhoxK4WEj49ed6XBYRAcnJ\n1u3YMTh+HPPYMXamHuS37KOs8jDYUL0G20PC2Vu7PnVSj9Bmz2467N3L5bv3cnniHuonHcLIzra2\nUTsVyry9rb8ZHx8rPHl4WGHqUl8NwwpqTufp29k/n33Lzz+9E8ipr3l5p0PbmeHtYt/7+p6u+dTt\n7GNFOeeMY04vL/oM+gdx61/lBP7k4A2YpZpscao3//LLL69aIc1utVdW2Xk5hP70If0D/fi8280l\nayQ31/rE/8ADcM89ri2wknp5/U/8fX8G9wY5mdLx+ks/oCzk5FhDmv/6F/Tta4Wz5s3dU0sZWpu8\nhduXf8Vqr+Z09D7BrJi+RFULdndZRXNqiHTlSmu9trg4q/etRYvToS062pq8c3KRX7tYkbSB1zf/\nzvz0TFL8mlLLPMrVQdV5IqoLbWuUcMJSRZGbawXu5OTTX5OTcR4+zMdvfoF5uDUhpBBKMqEkU9ux\nlwCHiRESAqGh1mzpGjWgWrVzbwEBUK0aqd6w0DzGIuMEcQ7Y7hXAkYBwHA4PauelEuVl0LlmLQbU\nb8kVoQ1KvuSSK5jmucHt1NfzHTuz5+/U7VTovNixIpyTezyD7KPH8XKaeJFPnmHg4e+Hx6mt5y4R\n+nK9PNkR4M2Waj5sC/Bhh78vuwP8+frR8ZUipHUBxpmmOeDkzxcc7hw7dmzBz7GxscTGxpZnqVXK\n/P0bGLB+G1+3jGRw/XYla2TdOmtx27g4XQx9CX/58xM+Pu7HpAY1ebiFG7Y6M0346ivrWrOoKHjl\nFWs5lUrujwNruWvVfLb6NKej13He79CbtjVqubus4svKgtWrT4e2hARr54SwMGvttjNvzZuX20QP\np+nkh8TlvLU9jl8zTTL9m9DQTOX6WrV5PKoL9f2qlUsdLpeXZ+0Je6qH66zgdd5jGRlQs6YVtkJD\nC932nchlyhcJrEvqyBGqU73RRl6c8gDtuncp9UhEdl42Sw5s4PuDW1iWdoStOU6OOIIxfUKpnpdG\npGce7apXp1NwOLG1m9CyWlDJR1AuwZUTGlw9OaKgPdMkunVrHCeDYW5mBpuO7GV96gG2ZJycqZvv\n5JDhwRFPf9J9a5DpH4xvdiYBy5fimbCawLxcgnLzWDZ/caUIaR7AZqyJAweA5cDNpmluPOs89aSV\ns1v//JT/pedwqPdwqnmV8H/qL754eksgDXueIy8/jx6/vM9KZwjftm1D3/CW5V/EqlXwt79ZbySv\nvQb9+pV/DW726/61jFr9Mxu9mtDe4yjvd7iKmOC6Ln+ecp11l58PO3ZYH5bOvG3fDuHh1vB1kyaF\nvzZubPXKlEJOfg6ztizkjY1xbPQKwenfkBZGOrfVjeTBJjHU8PK+dCNnKPPfWX6+tWfrxQLW2ceP\nHoXg4MJh6zzhq9CxGjUuOuO9PP82TNNkY8pOvtm3nl9T9rMpM5skpxcZnjXBOxj//HRqGzk08fWm\nffVg2gfV4fKQ+jTxD8S71EuDxAKl21PUVW3lOPNZl7afhJR9bDx2mG0ZR9mbnc2hvHxSnV5kOgLI\n96qBkX8C3/xjBBo51PaA+j7eNPEPpFVgKNHB9WgXXBf/8yyNVGn27jy5BMcbnF6C46XznKOQVs6c\nTid1588g0svJ0t4lHLLMy7NmAd55Jzz4oEvrq+jSs4/RYcFUkr3qsLzzVbSsUTar8l9QUpI1EeD7\n760FiO+6q8rNIjzbykNbuC9+PgmOCFo5UpnYqit9w5u5pG3bbHydm2vtn7ptmxXYTn3dvt0KddWr\nQ716ULeudTv1fb16VrgLDbUCSrVqBR+80k+kM3nTAj7cv4uNZgimRw2MdUl4/OxLy8PxzPjvfWX/\nZpyfb010SUkp2u1UCEtPh6CgS4esM29BQZVyiaF8Zz5bUnexKGkby1L3s+H4MfbkmqThywnPGuAd\ngpczk0Azi1oeThp4exHhF0AD/+pE+gfTrHoojQKCCfXyKtQb58plLi7W1sqVE8k0TXZnpLL5WBLb\nj6WwKzONfScyScrJJjnPSbrTIANvsh1+5HtUw5F3DF/ncQLJppYH1PfxoZF/dVpUq0m7oHCiQxoS\nWMJOikoT0opCIc091qXup/2KZYzKP8Rb/e8t2ae7zZutdaB++cUabhF2pO0hZskX+PrXYU33IdTy\nKc9hSFoAABdeSURBVMcZnPn51h6aY8da4XnMmEo3W7C01qUkcv+q7/nTGUYt8yhPNozk8ebdS9y7\nYfeV4At6cZxOouvVw3HwoLU49b59sH//6e8PHCgIObl5eSzsFMOcKzrz3eWd8cKk/7bNNJ63lOp/\nNiLb9CcLP7LwJaTuD4z++wgcHh5WsDvfLT//9PVIOTk4s7OZOnkeqUmx+JJNNY5TjWPUq76cK9o3\nxDh5ET3Hj1sX1J84YfVW1ax5/ltwcOHvTwWu4OD/b+/Ow6OsDj2Of89kY5KQAJKwyx4KAULCYimI\nRcVSxGq1QNXWVn24uFZ7ey22ttUrdUGtrU8XV+JaLSqWzargtWytyBZWwaAssiUs2cg6mZlz/3hH\nRVlKwkzmneT3eZ55CC/JOyeHk8xvzupMipcT+qxtBGyAjN4dWV+2h/WlB9hWWcqumhoO+oNUWA/V\nJOLzJBNMSIP41sQFakiytSThx9T7KC3yYWu8UOeBGg/UxBFnizh3dDJntW0bOuGB0CkPznPXBgLU\nBAPUBgPUBYPUBS1l1VV8stuPTUoFrweS4iApERKTIMELQR/GX0VisAqvrSPNBGkbB5mJCXRK8nK2\ntzW9UtvSNzWDnHbdaB3B1fMKaRJRn72L3dxrAP7r25L1yGL+9rufNu6df36+s1Jw1aoWfwj18r3r\nGFfwPn3TOrJq9GV4m/IFYvVqZzFHcrKzQECh+ZSO1FZw+7oFvFph8Zg4prSJ5/dDJtA2qWFDgm4+\nU/F0e6ustawu3swfC//FO2WVHPb25axgHROsn1uSMhhWWcvOjRt5+N79eHwD8IYimpca0uI3MuXS\nVmS0b+/MfzzRIz7+i5V8iYnsP3KEp5+voqp+CD4SOUprKkmlPmkjD/yxB/2GDXN6/VJTnT+93mbZ\nuxVNjen9rfXXUlR5kO0VxeyoKqHEV8O2Pbv569/LCCR1BW89eAPgDWISSunTN5Ekr/fz81KDFizO\nI9EYkjyGJOMhyeOhVZyH2oqjLFtcTaCiL5R7oSwVytJIqlnLvGf7csHXzyPe444RgcaENKy1Mflw\nii5NJRAI2CFDbrUQcH6DPjDH8uwjdlDejTYQCDT8hsGgtZMnW3vzzeEvbAx5bst8G7cw31608h/W\nHww23ROXlFh7443Wduhg7fPPO/8fctr8Ab+9f9M7NuOtZ6xZNNcOXvSUnb1z1Wl//Zo1a2xy8pzj\nkkly8ut2zZo1ESz5qR33c4614FwLBAK2zl9nX/losT3/nUetd94frHn3H7bru7PtDQWLbWFl+XH3\nC+f36dY6ayn+U9toDveKtFBuaVDW0dsMOS0FBQWhd0+hJnPPd6B1DzZPSP58cmuDGOMMsS1cCPPn\nh7OoMcFay13//jPX76nmum5ZvD1iPHFNtZDijTeclZrWwtatcM01WsTRQHGeOH4x8CIOjr+epYMH\nkJ6QyFWFe0l56zkuX/48BYd3nvLrc3NzycpaAgSPuRokK2spubm5ESz5qR33cw6QdoDN3dPIemMG\nyW89zw/2+Dic2JV7vjaKI2MuZM8Fk3l8yIX0TTl+iDyc36db66ylOGHbwENh4XkNfg3weDzk509j\nyJDbSU6eQ3LyHHJybiM/f1qDh/rDeS83ckcfoMSeuniYPgb75wSe2ruYJxszPNOmDbz8Mlx+OeTl\nQdeu4S+nC/mDfq58+9fMTfg69/buxS/7DGqaJy4qgltucU6AePVVZ16gnLFzO/RlWYe++IIBHt66\nlPz9nzJ0/Yek+RYxIa0VP836OsMzv7y33GcvLNddd/uXDr7Oz78h6i8sNqEGhv8DvrkNBgahQ1/8\nNTm0a1XLg4O/waUde572flrh/D7dXGfScLm52axd+4djVrA+1uj/x3Dey200J01Oy8kmOne6ZQYH\nJ/Tmrf7dGdfj3Mbd/KGHYM4cWLYs5nZJb6jy2nIuWHAHm9pdyqz+A/lB5+6Rf1Jr4cUXnT3PrrsO\nfvObFj8PMNJKfbXM3LacV4qL2OM5i6Ta/eQl1nF1l95c1WskbVqlA+7YJ2pP+R7m7l7FwqLdFNTU\ncyi+E8S1h+3VsLILvDeQIR3vOqMFDW74PuXMuH2xSyzQwgGJqGPPIAPnXeyzz97Az+uLWPHpCrZd\n+GO6t2lE6LAWJk1yVlk99VSYS+0eu8p2ce7CX3Oky9UsyBnGBe2aYEf1Tz+FadOcVXizZjm7z0uT\nqgkEyN+5jpf2fcwGXyI1JJBeu5vBSZYL2ndiYpcB5GT0O6PJzaczoTsQDLCjdAdLDmxh2eG9rK88\nyif+BGpT+pDkiSMr3sd5bduTW+PlsRtfYvtHX/45b/KtQcR1TvYaoLZxehTSJOJO9C62KhCg94pF\nxO+fx7YrHiE1sRG7hh89CuecAz/9KUydGuZSR9/KvSv51ntPYM6+iiVDRzIk0sfzBIPOnL/f/AZu\nuw2mT3fOupOo21FVwYu7N/DWof185DOUx6VDbTFt/IfpGh+kZ6skslPbktu2E31aZ9AhpT1nec8i\nKf7EvczBYJC8obex4cP7ILkEkg9Dh2Iyh/+VkZOGs7uuln31cMSThk3pSSJBOplavuZtxUWZ3bi8\ncxbdv9Kzqt4qORm1jcZTSJOo2VZVRe4HyxhRtoB/fvdPeEwj908791xYsMAJbM3E7M2vcu3GFbTr\ndhkrho2kR6SHGrdvd85H9fmc3rMBAyL7fHJGfMEga8sP84/9hWw4WsLO2loO+A0VJpl6TytMoBrr\nK8X4q/AQwGMDxBF0tidwrhAkHuJbQ0IbiPeCrw4qyumRFqBfejoDWrflgsyenNMmg/aJDdvhX0TC\nQyFNomruwSKmbPyAHwVW8tS3HmjcTebPh5tvhvffj/mFBNZafrv8AWYe8tOzy1jeyzuHjEi+QPr9\nzt5zM2fCr34Ft96qTTljXMBaSuvrOejzUVRXRZXf5zwCPjx4SIlPZP/O3Uy/rQTfwW9BRQIcjYeg\nxxV7ronIFxTSJOp+9fE2Hi18n7vTS5k+6r8bd5OHH3Ymui9f7uwYHoPq/HVcv/AmFsTnMbTzcBbk\nDCUlkoFp0yZnUUBaGjz9tHPmorQImtAtEhsaE9L00ythNaN3P8Z2HsiMokr+uvHlxt3kf/7HGfa8\n4gpnyC7GFFUWMfrFCbyVciETe47hnSHDIhfQfD7nOKfzz3cWCLz7rgJaC9Pc94kSacnUkyZhV+n3\nM3zN++z5+CX+PvJ7jOs9ruE3CQSc/dPatIHnnouZzVYLDhQw8Y3r8GX/lmu79WNmr96YSJV91Sqn\n96xXL3j8cefAa2mxNKFbxN003CmucaCujrzV71NV+GfevOAnnNu9EXuoVVfD2LFOL9H997s+qL22\n5TWmLv0dcQPv49c9s7i9W7fIPFFlpXMI+iuvOHPQvv9919eNiEhLp+FOcY1OSUn8M3c48X1u4pK3\n72X57uXHfU4wGGTt2rWsXbuWYDB4/E2Sk51joxYuhHvuiXyhGylog9yz5B5uXPU3zKCZ5A/IiVxA\ne/tt5xD0khLYvBmuvFIBTUSkmVJPmkTUv8rLmbhhHWy6k/kTH/m8R+10Nt/83MGDTm/apEnO/CsX\nKast45q5P2J9wgACnS5m4aAcciOxB9rhw84ecitWOPufXXRR+J9DREQiRsOd4kpzDx3i+q2bsRvv\nYO4ljzK62+iGr0Y7eNAZ+pwyxdmg1QXWHVjH916/kqT+d5HSpj/zBw2mc7iPtbLWOd/0Zz+Dq6+G\ne++FlJTwPoeIiERcY0KaDliXiLssI4N6m80NPMxlC37G9O7fC/WgHRvGPBQWnkdBQcGJ93XKzIT3\n3nOCWnW1M0ctShOjrbU8s+4Z7lz2IJnDHye7TWde6N+f5HCv4Ny+3dnr7MABZ4Pf4cPDe38REXE1\nzUmTJjEpM5MnvpaNJ+cRHto1h/rh84AG9oR26OAcwr58udOrVFcXkbKeypHqI0x+fTIPblpA4vBn\n+V6XvryanR3egFZVBXfdBSNHOsO8a9YooImItEAKadJkJmVm8ni//nhyf4cZ+yFMuAVMIPSvQbKy\nlpKbm3vqm7Rv7+wF5vc787JKSiJS1hMtalj0ySIGP5HD4bZjqOw7nWf6ZzOjZ0884Zq4by289hr0\n7w+7dsGGDfDzn+vMTRGRFkrDndKkJmVmkmAM1/ofJi7uUeo65pC48Db6ddpIfv4Np7e3k9cLs2c7\nh4aPGgVvvOEEmzD56qKG3gOeJvvWoywrXU3W6Bco86TxfnY2vcJ5BueWLc5B6AcPOqctnHde+O4t\nIiIxSQsHJCo+qKjgu5s30evQB2zbej8vXP4CE7ImNPxG+flOWJs5E6699oy3ozjuiJ1+82DCraQm\nDaHtJdOZ2D6D3/XujTdcw5v79jkrVufPd87bvOkmiNd7JxGR5karOyWm7Kqp4eJNm+gXV82qZT/i\nqoFTmDF2BknxDVwhuWWLs+ozJ8fZeT8trdFlWrt2LWPG7KY6YRh8+yeQuR2S/oQ5x/L7rincds45\njb73l5SWOmeUPvkkTJ0Kd97pnK4gIiLNkjazlZjSw+vl33l51CS0o+u5r7GprIicJ3JYumtpw26U\nne0ckZSa6gS1hQsbXaYyXxm+sc/BtFxI+CYMyYf09rS6tYzR4ejhOnLE6THr0weKi2H9enjwQQU0\nERE5jkKaRFV6fDxvDhrEpRkdKOg8jYmjHuGHf/8hk1+bzI7SHad/o+Rkp1fqySedPcUuuQR2nP7X\nl9eWM2PpDKYsn0Ja58OQsAguGQEv9IRfDaBfhyX/eVHDqezZ4wzLZmU5885Wr4ZZsyBSJxOIiEjM\nU0iTqPMYwy+6d2feoEG8WZ/J0AsX0D1zGMOfHs71866n8Ejhad0nGAyy9qyzWPfccwS/8Q0YMcIZ\nRiwqOunXbD20lZvfvJkej/XgwyMfc8Nli7FXP0BGwjq8NxaRvGo5OTm3k58/reEHVlvrbBcyaZLT\nw1dbCwUF8NRTzqHoIiIip6A5aeIqtYEA9336KY/v28eNHdvj3zObZ9b8ibxOeUzNm8rFfS/Gm3D8\nqsoTHTP14v3fYeCbc50d+ydNgjvugD59KKosYsFHC5i9ZTabD27muqHTSD17Cn8pLmNo69bM6NmT\ngcnJFBQUAJCbm9uwgPbpp87h5y+9BD6fsyHtNdec0Vw5ERGJbVo4IM3Gx9XV/HLnTpaXl3N7l06k\nl65gzqYXWL1vNef3PJ9xvcYxsttIBmQMIN7En/KYqcq9Oyh5+H/JfPENNnRLZNZAH4GLv82IQVM4\nlJbHM0UHGZSSwj09ejC8MUHKWti2DRYvhjlznIPPr7gCrroKxoyJ2skIIiLiHgpp0uysP3qUmXv2\n8E5JCZMzMpiYnkRp8Qr+b+e7rN6/mh2lO8hIzGDf1nSC1X2hPhk8fkiowtN6G217llIbrGVwh8GM\naT+My3eksm/HUZ7r2pXlgwczqbyc/+rYkaH9+ztHT53OFh4VFU4Q27jRWbCweLETxC66CCZOhPHj\nIdxneIqISExTSJNma29tLS8UF/NScTGlfj/j27VjVFoaA72JfLJhKdfftJE6BkNCDQTjoT6ZRP8W\nnnl2CB0GjGBNZSXLy8v5V3k5Q1u35gfp6UzZsIHUpUth5UooLHROMcjKgi5dnA1zW7VywlZ1NRw+\n7DyKi50/s7Nh8GDIzYVx46Bv3zPeo01ERJovhTRpET6pqeGdkhJWVlSw9uhRdtTWUl9RQ6AsHeri\nIM6CN4BpV0375FZkp6SQl5rKqPR0xrZpQ9uTHbN05IhzqPn+/c65oLW1ziM52TmOKiPDeZx9NoT7\nMHUREWnWFNKkRQpayz/XbebWO19m54FzIAA9O73PrEd+yMi8gdEunoiIiEKatGzBYLDxKzJFREQi\nSCFNRERExIV0LJSIiIhIM6GQJiIiIuJCCmkiIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIiIuJCCmki\nIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIiIuJC\nCmkiIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIiIuJCCmkiIiIiLqSQJiIi\nIuJCCmkiIiIiLqSQJiIiIuJCCmkiIiIiLuSqkGaMecgYs9UYs94YM8cYkxbtMomIiIhEg6tCGrAI\nyLbWDgG2A7+IcnnkGEuWLIl2EVoc1XnTU503PdV501OdxwZXhTRr7bvW2mDoryuBrtEsj3yZfqib\nnuq86anOm57qvOmpzmODq0LaV1wHvBXtQoiIiIhEQ3xTP6ExZjHQ4dhLgAXustYuCH3OXUC9tfbl\npi6fiIiIiBsYa220y/AlxpgfA1OB8621daf4PHcVXEREROQUrLWmIZ/f5D1pp2KMGQ/cAYw5VUCD\nhn+jIiIiIrHEVT1pxpjtQCJwJHRppbX2pigWSURERCQqXBXSRERERMTh5tWdJ2SMGW+M2WaMKTTG\nTI92eVoKY8wuY8wGY0yBMWZVtMvTHBljZhljio0xG4+51tYYs8gY85Ex5h1jTHo0y9jcnKTO7zbG\n7DXGrAs9xkezjM2JMaarMeY9Y8wWY8wmY8xPQtfVziPkBHV+a+i62nmEGGOSjDEfhF4vNxlj7g5d\nb3A7j6meNGOMBygELgD2A6uB71trt0W1YC2AMWYHMNRaWxrtsjRXxpjRQCXwgrV2cOjaTOCItfah\n0JuSttbaO6NZzubkJHV+N3DUWvtoVAvXDBljOgIdrbXrjTGpwFrgUuBa1M4j4hR1PgW184gxxiRb\na6uNMXHAv4CfAFfQwHYeaz1pI4Dt1trd1tp64G84jU0izxB77SWmWGtXAF8NwZcCz4c+fh64rEkL\n1cydpM7Bae8SZtbaImvt+tDHlcBWnE3L1c4j5CR13iX0z2rnEWKtrQ59mISzSNPSiHYeay+6XYA9\nx/x9L180NoksCyw2xqw2xkyNdmFakExrbTE4v2yBzCiXp6W4JXSG8DMaeosMY0wPYAjO6TId1M4j\n75g6/yB0Se08QowxHmNMAVAELLbWrqYR7TzWQppEzyhrbR4wAbg5NEwkTS925ifErr8AvUJnCBcB\nGg4Ks9Cw2+vAbaHena+2a7XzMDtBnaudR5C1NmitzcXpKR5hjMmmEe081kLaPuDsY/7eNXRNIsxa\neyD05yHg7zhDzxJ5xcaYDvD53JKDUS5Ps2etPWS/mKz7NDA8muVpbowx8Thh4UVr7bzQZbXzCDpR\nnaudNw1rbQWwBBhPI9p5rIW01UAfY0x3Y0wi8H1gfpTL1OwZY5JD78IwxqQAFwGbo1uqZsvw5Xki\n84Efhz7+ETDvq18gZ+xLdR765fmZy1FbD7d84ENr7WPHXFM7j6zj6lztPHKMMe0/Gz42xniBcThz\nARvczmNqdSd8firBYzgBc5a19sEoF6nZM8b0xOk9szgTIP+qeg8/Y8zLwDeBs4Bi4G5gLvAa0A3Y\nDUy21pZFq4zNzUnqfCzOvJ0gsAuY9tk8EjkzxphRwDJgE87vEwv8ElgFvIraedidos6vQu08Iowx\ng3AWBnhCj9nW2vuMMe1oYDuPuZAmIiIi0hLE2nCniIiISIugkCYiIiLiQgppIiIiIi6kkCYiIiLi\nQgppIiIiIi6kkCYiIiLiQgppIiIiIi6kkCYiIiLiQgppIiIhxhiPMWarMaZztMsiIqKQJiLyhaFA\nO2vt/mgXREREIU1E5AtjgfeiXQgREdDZnSIiGGMuA84DrgRWA9uBJ6y1hVEtmIi0aAppIiKAMSYR\nKAFyrbXbo10eERENd4qIOEYB5QpoIuIWCmkiIo4LgaXRLoSIyGcU0kREHBcCSwCMMaNDw58iIlGj\nkCYi4hgIfBAKZ9+w1vqiXSARadm0cEBEBDDGPAzUA4eAp6y1VVEukoi0cAppIiIiIi6k4U4RERER\nF1JIExEREXEhhTQRERERF1JIExEREXEhhTQRERERF1JIExEREXEhhTQRERERF1JIExEREXEhhTQR\nERERF/p/BgabFaz9QPEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106fd7e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t_train, y_train, 'o', label='data')\n",
    "plt.plot(t_test, y_test, label='true')\n",
    "plt.plot(t_test, y_lsq, label='lsq')\n",
    "plt.plot(t_test, y_robust, label='robust lsq')\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$y$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We clearly see that standard least squares react significantly on outliers and give very biased solution, whereas robust least squares (with soft l1-loss) give the solution very close to the true model."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
